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This  technical  report  consists  of  three  parts.  The  central  problem  is 
the  extrapolation  of  band-limited  signals. 

In  part  I,  several  existing  algorithms  for  band-limited  extrapolation  are 
compared:  Two-step  procedures  appeared  to  give  better  reconstructions  and  require 
less  computing  time  than  iterative  algorithms. 

In  part  II,  five  basic  procedures  for  iterative  restoration  are  unified 
using  a  Hilbert  Space  approach.  In  particular,  all  known  iterative  algorithms 


for  extrapolation  of  band-limited  signals  are  shown  to  be  special  cases  of 

7/v?- 

Bialy's  iteration.  ~We~also  -ebtai^  faster  algorithms  than  that  of  Papoulis- 

O  irv  ■*  at 


Gerchberg. 


In  part  III,  the  extrapolation  problem  is  presented  in  a  more  general 

f'fQ.SSn  /? 

setting:  Continuation  of  certain  analytic  functions.  Me  present-  two-eetrpg"  ^ 


procedures  for  finding  the  continuation  of  these  functions.  Some  new 
procedures  for  band-limited  continuation  are  also  discussed  as  well  as  the 
case  in  which  the  signal  is  contaminated  with  noise. 
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Numerical  Comparison  of  Several  Algorithms  for 
Band-Limited  Signal  Extrapolation 

Thomas  S.  Huang 
Jorge  L.C.  Sanz 
Hong  Fan 
Jamal  Shafii 
Bin-ming  Tsai 

Coordinated  Science  Laboratory 
1101  West  Springfield  Avenue 
University  of  Illinois  at  Urbana-Champaign 
Urbana.  Illinois  61801 

AP.SffACT 

We  present  some  computer  simulation  results  on  the  band-limited  signal 
extrapolation  problem.  First,  the  performance  of  several  existing  algorithms 
are  compared  for  the  noise-free  case.  Then  ve  describe  some  modifications  of 
these  algorithms  for  computing  the  extrapolation  when  the  given  signal  is 
contaminated  with  noise.  Computer  simulation  results  for  both  the  noiseless 

and  noisy  cases  are  included.  From  these  results,  the  following  preliminary 
conclusion  could  be  drawn:  Two-step  algorithms  appeared  to  give  better 


reconstructions  and  require  less  computing  time  than  iterative  algorithms. 
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I.  BflBflBggriflH 

The  band-limited  signal  extrapolation  problem  was  addressed  by  several 
authors  ([1]  -  [10],  among  others).  Soae  algorithas  for  coaputing  the 
extrapolation  were  also  given.  But  probably  the  aost  well-known  algorithas  in 
the  engineering  literature  are  those  of  [2],  [3]  and  [6],  [8].  Although  soae 
nuaerical  exaaples  were  given  in  [6],  a  nuaerical  comparison  between  both 
algorithas  seeas  not  to  be  available  in  the  easily  accessible  literature.  In 
this  paper,  the  nuaerical  performance  of  several  existing  algorithas  are 
coapared  by  means  of  computer  simulation  exaaples  (Section  II).  Then  soae 
aodif icatious  of  these  algorithms  are  proposed  tot  getting  the  extrapolation 
when  the  given  signal  is  contaminated  with  noise. 

Let  ns  recall  what  is  meant  by  band-limited  signal  extrapolation.  Assume 
that  g:  B  ->  C  is  an  fl-b and- limited  finite-energy  signal,  i.e. 

g(»)  «  0.  u  f  [-fl.fi], 

where  g  denotes  the  Fourier  transform  of  g.  If  we  are  given  a  piece  of  g,  g: 
[-A,A]  ->  C,  we  will  be  able  to  recover  g(x)  for  x^  [-A.A],  because  g  is  an 
analytic  function.  Band-limited  signal  extrapolation  is  the  problem  of 
computing  g(x)  for  all  x  from  the  known  values  g(x),  x€[-A,A] .  Several  basic 
models  relevant  to  the  solution  of  the  extrapolation  problem  were  given  in  [1] 
(also  [9]). 

Probably  the  most  well-known  technique  for  solving  the  problem  is  given 
by  the  iterative  procedure  ([2], [3]): 

80  ■  initial  O-band-limited  approximation 


gu  "  *incQ  *  <J[-A,A]«  +  <I"J[-A,A]>  *u-l> 


(1) 


where  J is  the  truncation  operator  to  [-A,A]  and  I  denotes  the  identity. 
Some  generalizations  of  this  procedure  were  given  in  [4].  The  numerical 
computation  of  (1)  can  be  accomplished  by  means  of  the  two  following 


techniques : 


(i)  Implementing  the  convolution  by  using  FFT 
(ii)  Sampling  the  iterative  equation 
b0  “  initial  guess 

hu  “  hu-l  +  *  J[-A.A](*  -  Vl>  <*> 

Note  that  recursion  (2)  is  equivalent  to  (1) . 

Technique  (i)  leads  us  to  the  following  discrete  recursion: 


N-l 


*u+1<j)  -  2N  Z  °u<*>  .  e2aija/N  ,  jet-N.N-1] 
m— N 


»— N,  k  «  1  |B|  i  N-l 


«.<■> 


(3a) 


(3b) 


N-l 


I  Ptt(j>e“2nijB/N.  I»l  <  kfl 


j— N 


Mj> 


f  «<J)> 


I j I  <  [rJ 


l.Tu(j),  j— *  <>*■  I*]  i  Ijl  i  N-l 


(3c) 


In  formulas  (3b)  and  (3c)  A  denotes  the  distance  between  consecutive  samples 

of  g  inside  [-A,A] .  It  is  to  be  noted  that  A  should  be  chosen  so  that  [41  • 

A 

.  The  convergence  of  this  procedure  was  shown  in  [5],  Some  relationships 


between  the  eolation  to  (3«)-(3b)-(3c)  end  the  eolation  to  the  extrapolation 
problem  were  given  in  [1].  In  Appendix  B,  we  will  ehow  that  the  Unit  of  the 
iterative  algorithm  (3a)-(3b)-(3c)  can  be  obtained  by  means  of  a  certain  two- 
step  procedure  ([6]).  (A  related  diecaeaion  is  given  in  [11]).  Naaerical 
comparison  between  both  algorithms  will  be  given  in  the  next  section. 

On  the  other  hand,  technique  (ii)  originates  the  following  iterative 
procedure: 

Su<J>  -  su-l<i>  +  *  l  sinc0[(j-k)A]*<g(kA)-S  (k))  (4) 

k€[-A,A]  “  1 

This  recursion  was  shown  to  be  convergent  in  [6].  It  was  also  shown  ([6]) 
that  the  limit  of  the  procedure  can  be  computed  by  aeans  of  a  two-step 
algoritha.  The  relationship  between  this  discrete  technique  and  the  solution 
to  the  extrapolation  problem  was  discussed  in  [7].  In  the  next  section,  the 
solution  given  by  this  two-step  procedure  will  be  compared  with  the  two 
techniques  mentioned  earlier. 

If  the  piece  of  the  signal  g  is  contaminated  with  noise,  the  situation 
will  be  completely  different.  First,  the  extrapolation  problem  will  not  have 
any  solution  (unless  the  noise  is  also  Q-band-limited) .  Therefore,  the 
problem  has  to  be  restated.  Several  attempts  were  made  in  this  direction 
(111. [6]). 

Even  though  equations  (1)  and  (2)  are  no  longer  convergent,  the  discrete 
iterative  procedures  (3a)-(3b)-(3c)  and  (4)  are  indeed  convergent.  However, 
the  extrapolations  obtained  by  means  of  the  two-step  algorithms  are  of  poor 
quality  in  the  presence  of  noise.  This  shows  that  some  sort  of  stopping  rule 
is  necessary  for  the  successful  application  of  the  iterative  procedures  in 


order  to  prevent  the  aoise  from  props gat log  too  aech  into  the  reeoaetraetioas. 

Te  also  present  soae  heuristic  resalts  obtained  by  modifying  the  two-step 
algorithms  to  cope  with  noise  in  the  given  part  of  the  signal.  This  will  be 
done  in  Section  III.  The  naaerical  resalts  obtained  will  be  compared  with  an 
iterative  Wienner-type  procedare.  The  convergence  of  this  iterative  algorithm 
is  proven  in  Appendix  A. 


In  practice,  the  case  in  which  g  is  not  contaminated  with  any  noise  is 
not  interesting.  However,  if  no  algorithm  can  perform  well  in  the  absence  of 
noise,  there  will  be  no  hope  to  solve  the  problem  in  cases  where  some  noise  is 
present. 


As  it  was  pointed  oat  in  the  Introdaction,  three  algorithms  will  be  ased 
for  the  extrapolation  of  g.  Another  algorithm  for  the  noisy  case  will  be 
presented  in  Section  III.  It  can  also  be  ased  for  the  noise-free  case;  this 
will  be  done  in  Section  III. 4.  In  oar  examples,  the  function  g  will  be  given 
by  the  formula 


The  Fourier  transform  of  g  is  shown  in  figure  0.  Note  that  g  is  band-limited 
to  [-1,1]. 

We  will  use  the  following  values  for  A;  Aj  .  i,  Jij  .  1/2  and  A3  =  1/4. 

It  is  worth  pointing  out  that  there  exists  a  fundamental  difference 
between  the  case  A^,  and  the  cases  and  A3 .  When  the  function  g  is  sampled 


over  z  *  (-1.1)  (say,  32  staples)  tad  DFT  of  these  staples  tre  used  for 
approxiatting  g,  the  two  peaks  in  the  frequency  space  are  still  distinguished 
(see  Fig.  la).  On  the  other  hand,  this  will  not  be  the  case  for  (-1/2, 1/2) 
and  (-1/4, 1/4) .  The  corresponding  plots  are  given  in  Figures  1(b)  and  1(c). 
All  the  extrapolations  will  be  coaputed  for  x  €  (-8,8). 

II. 1 

In  this  section  we  will  apply  the  discrete  iterative  procedure  given  by 
equations  (3a)-(3b)-(3c) .  For  all  the  cases  we  will  use  32  samples  in  the 
known  range  of  the  signal  g.  Therefore,  32  unknown  frequency  values  are  to  be 
sought.  For  A  «  1  the  length  of  the  DFT  will  be  256,  for  A  =  1/2  it  will  be 
512  and  for  A  »  1/4  it  will  be  1024.  Figures  2(a),  2(b)  and  2(c)  depict  the 
resolution  of  the  Fourier  transfora  obtained  after  10  iterations  for  (-1,1), 
(-1/2, 1/2)  and  (-1/4, 1/4),  respectively.  Figures  3(a),  3(b)  and  3(c)  show  the 
result  obtained  after  100  iterations.  Several  conclusions  can  be  drawn  from 
these  examples.  For  the  case  A  -  1,  some  negative  values  of  the  spectrum  have 
been  removed  and  the  two  peaks  are  also  more  clearly  distinguished  after  10 
iterations.  (Compare  Figures  1(a)  and  2(a)).  However,  the  improvement 
obtained  after  100  iterations  is  not  significant  (see  Figure  3(a)).  For  the 
case  A  **  1/2  and  A  *  1/4,  the  peaks  are  not  distinguishable  after  ten 
iterations  (see  Figures  2(b)  and  2(c)).  This  situation  is  the  same  for  A  = 
1/4  when  the  number  of  iterations  is  100  (see  Figure  3(b))  and  1,000  (figure 
4(b)).  On  the  contrary,  the  situation  improves  for  A  =  1/2  when  the  number  of 
iterations  is  increased  up  to  1,000  (see  Figure  4(a)). 
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I.« 


As  it  was  pointed  oat  in  Section  I>  the  limit  of  the  iterative  procedure 
(3a)-(3b)-3(c)  can  also  be  obtained  by  means  of  a  two-step  procedure.  This 
two-step  procedure  vas  also  given  in  [6].  However,  in  [5],  the  relationship 
between  this  procedure  and  the  algorithm  (3a)-(3b)-(3c)  was  not  discussed. 
Let  L  be  the  following  (2n+l)  x  (2n+l)  matrix, 

L(k.h)  e2nij(k-h)/M#  _n  <  k#h  i  n  (5) 

j—  n 

where  M  *  2N  +  1 .  The  matrix  L  is  positive  definite  and  therefore,  we  can 
always  compute  the  solution  x:  *xh>-a<i,£n  of  the  *7***“  of  equations: 

Lx  “  y  (6) 

The  two-step  procedure  consists  of  first  solving  equation  (6)  when  y  «*  (g(kA), 
k  =  -n,...,n),  and  then  computing  the  extrapolation  as  follows: 

zk*  *  M  ,  ^  2sij  (7) 

h— n  j— n  h 

where  -•  <  k  <  +  •.  The  extrapolation  is  the  limit  of  the  iterative 

algorithm  given  by  (3a)-(3b)-(3c) .  The  proof  of  this  fact  is  relegated  to 
Appendix  B.  This  way  of  proving  the  convergence  of  (3a)-(3b)-(3c)  is  simpler 
than  that  of  [3]  where  non-expansive  properties  of  certain  operators  were 
used.  Ve  have  chosen  an  odd  number  of  points  so  that  L  is  real.  In  our 
examples  we  will  use  33  points.  Figure  3  depicts  the  Fourier  transform  of  the 
extrapolations  obtained  by  using  this  two-step  procedure.  Figure  3(b)  shows  a 
much  better  resolution  for  A  m  1/2  than  that  obtained  in  Section  II. 1  (compare 
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Figure  5(b)  with  Figures  2(b),  3(b)  snd  4(a)).  Figure  5(c)  shows  the  result 
obtained  for  A  *  1/4.  Tbe  quality  of  tbe  reconstruction  is  good,  although 

soae  artifacts  are  present  in  the  lower  part  of  the  peaks.  This  exaaple  also 
shows  how  slow  the  iterative  technique  given  in  Section  II. 1  is  (Figure  4(b)) 
since  it  should  converge  to  the  saae  solution.  Ve  recall  that  the  result  of 
Figure  4(b)  was  obtained  after  1,000  iterations. 

II. 3 

Ve  now  consider  another  two-step  procedure  ([6],  [8]).  As  we  have  pointed 
out  in  Section  I,  this  two-step  procedure  will  give  us  the  limit  of  the 
iterative  approach  (4) .  This  iterative  procedure  (4)  and  some  of  its 
generalizations  were  found  to  be  very  slow  ([7])  and  no  numerical  examples 
using  them  will  be  included  here. 

Let  us  consider  the  matrix  K  €  E(2n+1) x(2n+l ) 


Ki.j,  .  a.iMflXjri).  , 


(8) 


Ve  now  solve  the  system  of  equations; 


Kx  *  y. 


(9) 


where  y  *  (g(kA),  k»-n,...,n).  If  x  =  (xh)-n<h^n  denotPS  the  solution 
(which  exists  and  is  unique),  the  extrapolation  is  computed  by  means  of  the 
formula 


n 

X 

h»-n 


sinnflA( i-h) 

it(i-h)  xh 


si,  -»  <  i  <  +® 


(10) 
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The  Fourier  transform  of  the  extrapolations  obtained  by  using  this  technique 
are  shown  in  Figures  6(a).  6(b)  and  6(c).  It  is  worth  pointing  oat  that  the 
resolutions  obtained  for  all  the  cases  are  very  good.  For  A  *  1/4*  the 
estimated  Fourier  transform  is  ewen  better  than  that  of  Figure  S(c)  since,  for 
this  algorithm,  no  artifacts  are  introduced  in  the  reconstructions  (see  Figure 
6(c)). 

To  conclude  this  section  we  would  like  to  point  out  some  differences 
between  the  two-step  procedures  (6) -(7)  and  (9) -(10).  Both  procedures  are 
intended  to  provide  approximations  to  the  solution  of  the  continuous 
extrapolation  problem.  However,  the  nature  of  the  two  approximations  are 
completely  different.  The  extrapolation  provided  by  the  two-step  procedure 
( 6) — (7 )  is  a  M-periodic  discrete  signal;  xfc>  — <k<+-  which  is  band  limited  to 
[-n.n] : 

N 

I  z  e-2nikj/M  .  0,  | j |  >  n,  M  =  2N+1 
k— N 

and  z^  m  g(kA)»  k  in.  On  the  other  hand,  the  extrapolation  given  by  (9)-(10) 
is  a  finite-energy  (and  therefore  non-periodic)  sequence  s^:  -«<k<+«,  which  is 
band-limited  to  [-QA.QA],  i.e. 

I  sv  e2nikw  -  0.  |w!  >  QA 

k—  * 

and  s^  -  g(kA),  1 k I  i  n.  It  is  known  that  if  A  ->  0,  then  the  extrapolation 
s]c  will  approach  the  continuous  extrapolation  g  ([7]).  On  the  contrary,  no 
similar  result  is  known  for  the  periodic  extrapolation  z^.  A  more  detailed 
discussion  about  four  basic  models  for  extrapolation  which  are  related  to  the 


two-step  procedures  can  be  found  in  [1] 


III.  THE  NOISY  CASE 

In  this  section  we  will  discuss  several  techniques  for  solving  the 
extrapolation  problem  when  the  given  signal  is  contaminated  with  some  additive 
noise,  the  noise  which  will  be  used  is  white  and  uniform  and  its  values  will 
range  on  (-0.005,  0.005),  on  (-0.05,  0.05),  and  on  (-0.5,  0.5).  They  will  be 
called  third  digit,  second  digit  and  first  digit  noise  respectively. 

In  what  follows,  we  will  compare  the  iterative  procedure  given  by 
formulas  (3a)-(3b)-3 (c)  applied  to  the  noisy  case  with  some  heuristic 
modifications  of  the  two-step  procedures  described  in  Sections  II. 2  and  II. 3. 
In  Section  III. 4  these  algorithms  will  be  compared  with  a  new  iterative 
procedure  designed  to  cope  with  noise. 

III.l 

We  have  mentioned  in  Section  I  that  the  iterative  procedure  3(a)-3(b)- 
3(c)  can  also  be  applied  to  cases  where  the  samples  g(kA),  -n  <  k  <;  n  are 
corrupted  with  noise. 

For  A**l,  the  algorithm  distinguishes  the  peaks  well  in  the  noise-free 
case  (see  figs.  1(a),  2(a)  and  3(a)).  Adding  second  digit  noise  to  the  known 
samples  causes  some  artifacts  in  the  reconstruction  (see  figures  7(a)  and 
7(b)).  This  is  also  the  case  when  first  digit  noise  is  used  (see  figure 
7(e))  . 

The  solutions  provided  by  the  algorithm  are  almost  the  same  as  in  the 
noise  free  case  when  third-digit  noise  is  added  to  the  known  samples.  This  is 
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•Iso  the  case  when  second-digit  noise  or  first  digit  noise  is  tdded  to  the 
original  signal  g  on  (-1/2, 1/2)  and  (-1/4, 1/4) .  Figures  7(c)  and  7(d)  shows 
the  resnlt  obtained  after  S00  iterations  for  (-1/2, 1/2)  and  (-1/4. 1/4) 
respectively  when  second  digit  noise  is  present  in  the  saaples.  Figure  7(f) 
and  7(g)  plot  the  corresponding  Fourier  transform  after  500  iterations  for 
first  digit  noise.  Coapare  the  results  so 'obtained  with  those  of  figures  3(b) 
and  3(c) . 

We  know  that  the  liait  of  this  procedure  coincides  with  the  extrapolation 
obtained  by  using  the  two-step  algoritha  described  in  Section  II. 2  (see 
Appendix  B) .  We  have  tried  the  two-step  procedure  given  by  formulas  (6)  and 
(7)  when  some  noise  is  added  to  the  saaples  g(kA),  “»  <  k  £  n.  The  results 
are  completely  wrong  due  to  the  presence  of  noise.  However,  the  iterative 
procedure  seeas  not  to  be  so  sensitive  to  the  noise.  This  is  due  to  the  very 
slow  rate  of  convergence  of  the  algoritha,  which  takes  1.000  iterations  to 
build  part  of  the  peaks  for  A  •  1/2  (figure  4(a))  and  probably  auch  aore  than 
that  for  A  =  1/4  (figure  4(b))  when  no  noise  is  present.  Therefore,  the  first 
thousand  iterations  will  not  be  enough  to  get  a  reasonable  approximation  to 
the  true  extrapolation  if  soae  noise  is  added  to  the  given  saaples.  It  is 
clear  that  any  attempt  to  speed  up  the  convergence  of  this  iterative  procedure 
will  also  propagate  the  noise  faster  than  the  present  algoritha. 

III. 2 

In  this  section,  we  modify  the  two-step  procedure  given  by  formulas  (6) 
and  (7)  to  cope  with  noise.  The  modification  consists  of  adding  a  positive 
number  y  to  the  main  diagonal  of  the  matrix  L,  obtaining  a  new  matrix  C,  then 
solving  the  system  of  equation  Ex  ■  y  and  using  formula  (7)  to  get  the 
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extrapolation. 

The  Motivation  for  this  modification  is  two-fold.  It  is  well  known  that 
L  is  an  ill-conditioned  matrix;  therefore  small  perturbations  introduced  in 

g(kA),  -n  <  k  i  n  may  produce  large  changes  in  the  solution  x^,  -n  £  h  £  n, 
and  this  may  overwhelm  the  extrapolation:  z^,  -N  £  h  £  N.  If  we  add  some 
small  positive  number  to  the  diagonal,  the  matrix  will  become  better 
conditioned  and,  therefore,  the  solution  of  the  system  of  equation  (()  will  be 
a  more  stable  problem. 

The  other  reason  for  such  modification  is  that  the  two-step  procedure  was 
shown  to  be  the  best  estimation  to  the  extrapolation  problem  when  the  noise  to 
signal  ratio  is  added  to  the  diagonal  of  the  matrix  K  given  by  formula  (8) 
(see  [6]).  It  is  obvious  that  K  and  L  are  different.  However,  some  related 
optimality  property  might  be  proven  when  L  is  used  instead  of  K. 

Figures  8(a),  8(b),  8(c)  and  8(d)  show  the  results  obtained  for  the  case 
A  ■  1  by  means  of  this  regularization  technique.  Third-digit  noise  has  been 
added  to  the  signal.  Even  though  the  spectra  are  different,  the  sensitivity 
of  the  reconstruction  with  respect  to  the  parameter  y  is  not  large. 

Figures  9(a),  9(b),  9(c),  9(d)  shows  the  results  obtained  for  A  *  1  when 
second-digit  noise  is  added  to  the  given  signal.  In  this  case,  the  sensitivity 
of  the  reconstruction  with  respect  to  k  is  more  evident. 

Figures  9(e),  9(f)  corresponds  to  the  first  digit  noise  case.  It  is  worth 
noting  that  the  case  given  in  9(e)  will  be  a  reasonable  reconstruction  if  the 
available  information  about  the  exact  band-width  is  used  before  plotting  the 


Fourier  transform. 


Figures  10(e),  10(b)  ead  10(e)  show  the  spectre  of  the  reconstructed 
sigael  for  A  ■  1  when  second-digit  noise  is  used.  The  result  obteiaed  is  not 
of  good  quality.  However,  soae  aore  e  priori  information  about  the  true 
spectrua  will  be  of  great  help.  For  instance,  a  positivity  constraint  applied 
to  figures  10(a)  or  10(b)  will  provide  a  auch  better  result.  This  is  also  the 
ease  when  first  digit  noise  is  used.  Figures  10(d),  10(e)  show  the 
reconstruction  for  this  case. 

Figures  11(a),  11(b)  and  11(e)  show  the  results  obtained  for  the  case  A  * 
1/4  when  second-digit  noise  is  used.  It  is  clear  that  the  sensitivity  to  the 
paraaeter  X  is  auch  aore  critical.  This  is  also  the  case  when  first  digit 
noise  is  added  to  the  signal  (see  fig.  11(d),  11(e)).  Once  aore,  a  positivity 
constraint  will  provide  a  good  reconstruction  if  X-0 .000003  is  used  (fig. 
11(a)).  Figure  12  depicts  the  results  obtained  when  third-digit  noise  is 
added.  Te  would  like  to  eaphasixe  that  if  positivity  of  the  Fourier  transform 
is  used  as  a  priori  information,  that  is  to  say,  we  set  to  zero  all  negative 
values  in  the  reconstructed  Fourier  transfora,  the  results  are  better  than 
those  obtained  when  no  extrapolation  is  performed.  For  instance,  fig.  11(f) 
depicts  the  Fourier  transformation  of  the  signal  plus  noise  in  (-1/4, 1/4)  when 
no  super-resolution  is  tried.  Compare  this  result  with  that  of  fig.  11(d) 
improved  when  positivity  information  is  incorporated  in  the  spectrua. 

III. 3 

Ve  have  pointed  out  in  III. 2  that  the  same  regularization  technique  can 
be  applied  to  the  two-step  procedure  given  by  formulas  (9)  and  (10)  (see 
Section  II. 3).  Some  motivations  for  such  techniques  can  be  found  in  (6] .  If 
we  denote  by  Z  the  matrix  E  +  XI,  where  I  is  the  identity,  the  technique 
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consists  of  solving  tie  system  of  equations  K  *  ■  y  (where  y  ■  (g(kA))  k  - 
-n.....  n)  and  nsing  equation  (10)  to  obtain  tie  extrapolation.  Figure  13(a) 
slows  tie  results  obtained  with  this  technique  for  the  case  A  «  1  and  when 
third-digit  noise  is  used.  Figures  13(b)  and  13(c)  show  the  results  obtained 
for  two  different  values  of  X.  where  second-digit  noise  is  added  to  the 
signal.  He  corresponding  plots  for  first-digit  noise  are  shown  in  fig. 
13(d),  13(e). 

Figures  14  and  15  show  the  reconstruction  obtained  when  third-digit, 
second-digit  and  first-digit  noise  is  used  for  the  case  A  -  1/2.  Figures  16 
and  17  present  the  corresponding  results  for  the  case  A  *  1/4.  For  this 
technique,  we  have  shown  the  plots  corresponding  to  the  best  results  obtained. 

Several  conclusions  can  be  drawn  from  these  examples.  For  A  -  1/4  and  A  « 
1/2,  when  first-digit  noise  is  used,  it  is  seen  (figs.  15(c),  15(d),  17(c), 
17 Cd)  and  17(e))  that  none  of  the  values  used  for  X  provide  a  reasonable 
result.  In  this  case,  the  main  conclusion  is  that  the  reconstructed  Fourier 
transforms  are  of  very  poor  quality.  (The  case  A  -  1/4  is  even  worse  than  A  * 
1/2)  For  A  ■  1  and  first-digit  noise,  the  results  look  much  more  encouraging. 
Taking  into  account  the  amount  of  noise  introduced  in  the  observation,  we 
conclude  that  the  reconstructions  are  acceptable.  As  we  have  remarked  above 
they  will  improve  if  some  positivity  constraint  is  used  to  remove  the  negative 
artifacts  (see  fig.  13).  Another  important  observation  is  that  the  results 
obtained  for  A  =  1/4,  A  “  1/2  are  much  better  when  the  technique  given  in 
III. 2  is  used  instead.  Comparing  fig.  10(e)  with  fig.  15(d)  and  fig.  11(d) 
with  fig.  17(c),  17(d),  17(e)  we  conclude  that  the  two-step  procedure  used  in 
III. 2  suits  this  numerical  simulation  better  than  that  of  section  III. 3  when 
noise  becomes  high  enough.  In  case  of  lower  noise  we  see  that  the  best  values 
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for  X  need  to  bo  smaller  as  tho  obsorratioa  part  of  tbo  signal  shrinks.  In 
addition,  tho  sensitivity  of  tho  reconstructed  Fourier  trass fora  in  torus  of  X 
increases  as  the  length  of  the  knova  part  decreases:  for  ezaaple,  figures 
17(a)  and  17(b)  show  that  the  result  may  change  drastically.  It  is  also  clear 
that  a  smaller  X  results  in  wore  positive  and  negative  overshoots  and*ringing. 
Again,  soae  a  priori  inforuation  nay  help  improve  the  reconstruction.  This 
will  be  the  case  if  a  positivity  contraint  is  used  because  the  negative 
artifacts  is  figures  13(d).  14(b).  14(c).  15(b).  15(c).  15(d).  16(a).  16(b) 

and  17(b)  will  then  be  removed. 

III. 4 

In  this  Section,  we  present  a  new  iterative  technique  for  the  noisy  case. 
Ve  will  sake  use  of  the  Periodic  Discrete  Prolate  Spheriodal  Sequences  (P- 
DPSS)  (see  Appendix  A) . 

Let  f(m)  be  a  periodic  bandlimited  signal  of  period  2N-1, 

f (k)  -  0  for  U<lk|<N-l  (11) 

where  f(k)  is  the  DFT  of  f(m),  and  U<N. 

Let 


y(m)  ■  f(m)  +  n(n),  for  0  <  B  <  D 


(12) 


y(m)  is  the  observed  sequence  of  length  D.  The  operator  T  is  defined  as 


T  x(m) 


x(n) ;  0  <  B  <  D 


(13) 


0  ;  otherwise 
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and  the  operator  B  is  defined  as  the  following 


B  x(b)  -  ID  FT  ^ 


x(b) ;  Ikl  <  D 


0  ;  otherwise 


(14) 


Also  for  convenience,  we  define 


(x(a)  ,  y(«))N  -  x(a)y(B) 

n-0 


(x(b)  ,  y(»))|j  »  2_  x(a)y(B> 
B*0 


(15) 


(16) 


For  the  given  noisy  observation  y(a) ,  we  can  expand  it  in  teras  of  0^(a) 
as  follows 


y(a) 


bk 


a€[0,D-l] 


(17) 


where 


bj.  »  1/Xj.  (y(a) ,  4k(a))jj 

*  1/kj.  (f(a),  dk(a))M  +  l/\k  (n(a),  4k  (b))d 

*  \  +  l/Xk-nk 


and 


*k  “  ( f (a) ,  0k(a) )fl 
nk  ®  (n(a) ,  6k(a) )D 

We  wish  to  find  a  set  of  Cj.'s  such  that 


*,(■) 


£  ckVk(.,  ;  0  1  ■  <  N 


(IS) 


ainiaizes  t he  aeaa  square  error  E  (f-f,,  f-f,)N. 

Usiog  the  orthogonality  properties- of  #k's»  we  eaa  write  the  aeaa  square 
error  as 

'(£<•*- ) 


ainiaization  of  the  aeaa  square  error  leads  to 


m  E  (akbk} 
k  E  (bk2) 


„2  2  _i 

+  a 
“k 


(19) 


where  we  have  assnaed  that  the  noise  and  the  signal  are  independent. 


Equation  (19)  is  siailar  ia  fora  to  the  faailiar  Veiner  filtering.  In 

practice,  a2,  a2,  and  X?  are  difficult  if  not  inpossible  to  obtain.  A 

*k  "k  k 

possible  approxination  is  to  assuae 


constant 


ck  -  1/(1  + 

kk 


(20) 


(21) 


Even  with  these  approxinations,  direct  solution  still  require  solving  for 

^k's.  Fortunately,  there  exists  an  iterative  algoritha  which  converges  to  the 
function  given  by  Equation  (18)  when  c.  is  given  by  (21)  . 
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Iterative  Algorithm: 

L*t  f0(«)  -  0 

Vl(")  “  Cl-n)£m(»)  +  BT  [BT[y(m)  -  fu(m)]]  (22) 

then  lia  f 
u->« 

The  proof  of  its  convergence  is  relegated  to  Appendix  A.  It  is  clear  that  if 
there  is  no  noise  present,  we  obtain  another  algorithm  for  the  noise-free  case 
by  setting  p  *  0  in  formula  (22) . 

Figures  18(a),  18(b)  and  18(c)  show  the  resnlt  obtained  with  p  *  0  after 
200  iterations.  In  fig.  18(a),  the  Fourier  transform  of  the  extrapolation  for 
A  *  1  is  shown.  The  result  can  be  judged  as  good.  However,  figures  17(b)  and 
17(c)  show  that  the  procedure  suffers  the  same  drawback  as  the  algorithm  given 
in  II. 1:  The  number  of  iteration  required  to  distinguish  the  peaks  is 
enormous . 

Figures  19(a)  and  19(b)  shows  the  performance  of  the  algorithm  for  (-1,1) 
when  second-digit  noise  is  added  to  the  given  signal.  The  value  used  for  p 
was  0.0.  Figures  20(a)  and  20(b)  shows  the  result  obtained  by  the  same 
procedure  when  p  *  0.01.  It  is  clear  that  the  new  value  chosen  for  p 
eliminates  the  artifacts  of  figure  19(b)  after  200  iterations  (figure  20(b)). 

Figure  21  show  the  results  obtained  by  using  the  procedure  for  A  *  1/2 
where  second-digit  noise  is  used.  Once  more, the  different  performance  between 
p  *  0  and  p  +  0  is  clear.  Figure  21  (b)  shows  that  the  artifacts  of  Figure  21 
8a)  have  been  eliminated.  However,  the  peaks  are  less  distinguishable.  The 
results  obtained  for  the  case  A  -  1/4  are  similar  to  those  of  section  III. 4 
and  they  will  not  be  included  here. 
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IV.  CONCLUSIONS  AND  DISCUSSIONS 

The  numerical  performance  of  several  algorithms  was  compared.  Two  basic 
assumptions  ware  aada:  tie  continuous  cut-off  frequency  of  tba  giwaa  signal 
is  known  exactly*  and  tba  Fonriar  transfora  of  tba  signal  was  supposed  to  be 
real.  In  tba  numerical'  axaaplas  glean  in  Sections  II  and  III  tba  iauginary 
part  of  tba  reconstructed  Fonriar  transfora  was  negligible.  It  turns  out  froa 
these  nnaarical  axaaplas  tbat  tba  non— iterative  techniques  (two-step 
procedures)  produce  better  results  than  those  provided  by  the  iterative 
procedures.  More  numerical  examples  are  in  order  to  verify  or  disprove  this 
conclusion.  Another  important  aspect  of  the  algorithms  which  has  to  be 
compared  is  the  number  of  operations  involved  in  the  procedures.  If  2n  is  the 
number  of  given  samples  of  g  in  (—A. A)  both  two— step  procedures  need  0(nA) 
operations.  This  is  because  the  matrix  involved  in  the  system  of  equations 
(6)  and  (9)  are  Toeplitz  and  therefore  0(n2) -algorithms  are  known  for  solving 
the  system.  On  the  other  hand,  at  every  step  of  the  iterative  procedures 
(3a)-(3b)-(3c)  and  (22)  at  least  one  FFT  is  needed.  The  length  of  this  FFT  is 
N«cn2  c  denotes  a  constant;  this  is  because  the  following  equation  has  to  be 
satisfied:  [^"Un.  Therefore,  we  will  have  0(n2log  n)  operations  per  step. 


satisfied:  Jj^J-n.  Therefore,  we  will  have  OUMog  n)  operations  per  step. 
This  analysis  shows  that  the  two-step  procedures  are  less  expensive  in  terms 
of  arithmetic  computation  time.  However,  the  performance  of  the  two-step 
procedures  for  the  noisy  case  depends  on  choosing  the  correct  parameter  A. 
Therefore,  the  relationship  between  the  optimal  A  and  the  noise  has  to  be 
further  studied. 

The  sensitivity  of  all  these  algorithms  to  changes  in  the  cut-off 
frequency  has  to  be  investigated  because  the  exact  value  of  the  highest 
frequency  might  not  be  available  in  practice. 
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Appendix  £ 

To  establish  notation,  sons  properties  of  Periodic  Discrete  Prolate 
Spherodial  Sequence  (P-DPSS)  [6]  are  listed  below. 

For  D  >  0  and  D  >  0,  we  can  find 


^o(a) «  ^(a)  /  •  .  .  «  »  ■  a  Z 

^0.  Iji  ....  such  that 


(i)  BT#^(a)  *  ;  0  £  «  <  N 

k  -  0.1,  .  .  .  ,  D-l 

(ii)  #k(«+N)  -  #k(«) 

(iii)  (^(a)  i  (a))jj  ■  8k  j 

(f*k(*)  ,  dj  (®) )  n  *  ly,  ij,  j 

(iv)  B  #k(a)  •  0k(m) ,  #k  is  bsndliaited 

(y)  #k(m) ,  k  ■  0,1,  ....  D-l  form  a  basis  in  the 

vector  space  of  sequences  of  length  D. 

B  and  T  are  defined  by  equation  (13)  and  (14)  respectively. 

Ve  will  now  prove  the  convergence  of  the  procedure  (22), 

Section  III. A. 


Let 


f  (m)  «  d.  0k(m),  0  i  a  <  N. 


(A.l ) 


by  substituting  Equation  (A.l)  and  Equation  (17)  into  Equation  (22),  we  have 
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r>  u 
► 


fc  I 


t  i 


dk,u+l  "  +  Xk(bk  -  dk>a) 

-  (1  -  <h*X2)) dk>a  +  x2bk 


(A. 2) 


so 


X?b, 


(l-(l-(xf  +  m))  ) 


k,u  ,2  'A'u"'Ak 
XJ+ji 


(A. 3) 


cbosing  p>0  such  that  1 1- (X^  +  |x)  |  <  1,  or  u  +  X2  <  2 

K  k 


we  have 


lim  f  („) 
u->« 


r  Xkbk 

t-  2  ^k(«n)  ,  which  proves  Equation  (23). 

kk+H 


a 

I 


L 


Agp?a4u  1 


We  prove  that  the  two-step  procedure  given  b  formal  as  (6)  and  (7)  in 
Section  II. 2  provides  the  same  extrapolation  as  the  iterative  technique  (3a)- 
(3b)-(3c) . 

Ve  know  that  the  matrix  L  given  by  (5)  is  symmetric  and  positive 
definite.  It  is  also  easy  to  prove  that  all  its  eigenvalues  are  less  than  1. 
Therefore,  the  solution  of  the  system  of  equations  (6)  can  be  computed  as 
follows: 


1°  -  0 

xu+l  -  Xa  +  (y  -  L  Xa)  .  a  >  0  (B.l) 

If  we  also  denote  by  L:  L(k.h)  the  matrix  defined  by  (5)  when  ke(-«»,+,“) , 
he(-N,N]  and  if  we  apply  L  to  both  sides  of  equation  (B.l).  we  will  obtain 

y°  -  0 

yu+1  -  yQ  +  LT(y  -  yu)  (B.2) 

where  T(y  -  yu) (m)  =0  if  m  £  [-n.nj  and  (y-ya) (m)  if  m  €[-n,n].  It  is  easy  to 
see  that  Lyu  *  yu,  for  all  u  >  o.  Therefore,  equation  (B.2)  can  be  written  as 

y°  -  0 

yu+1  -  L(yu  +  Ty  -  Tytt)  -  L(Ty  +  (I-T)ytt)  (B.3) 

where  (I-T)yu(m)  *  0  if  m  6[-n,n]  and  yu(m)  if  m  6[-n,n]. 

It  is  now  a  simple  exercise  to  verify  the  equivalence  of  the  procedures  (B.3) 


t'lguru  l.c.  DFT  of  31  if  Im  of  |(i)  In 
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ABSTRACT 

This  paper  deals  with  iterative  solutions  of  the  linear  signal  restora¬ 
tion  problem:  g=Af.  First,  several  existing  techniques  for  solving  this 

problem  with  different  underlying  models  are  unified.  Specifically,  the  fol¬ 
lowing  are  shown  to  be  special  cases  of  a  general  iterative  procedure  (Bialy 
1959  11])  for  solving  linear  operator  equations  in  Hilbert  spaces:  1.  A  Van 
Citter-type  algorithm  for  deconvolution  of  discrete  and  continuous  signals. 
2.  An  iterative  procedure  for  regularization  when  g  is  contaminated  with 
noise.  3.  Papoul i s-Gerchberg ' s  algorithm  for  extrapolation  of  continuous 
signals  ([2], [3]).  4.  An  iterative  algorithm  for  discrete  extrapolation  of 
band-limited  infinite-extent  discrete  signals  (and  the  minimum  norm  property 
of  the  extrapolation  obtained  by  the  iteration  [4])  and  5.  A  certain  iterative 
procedure  for  extrapolation  of  band-limited  periodic  discrete  signals  [5]. 
The  Bialy  algorithm  also  generalizes  the  Papoul i s-Gerchberg  iteration  to  cases 
where  the  ideal  low-pass  operator  is  replaced  by  some  other  operators. 

In  the  second  part  of  the  paper,  a  suitable  modification  of  this  general 
iteration  is  shown.  This  technique  leads  us  to  new  iterative  algorithms  for 
band-limited  signal  extrapolation.  In  numerical  simulations  some  of  these 
algorithms  provide  a  fast  reconstruction  of  the  sought  signal. 
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I.  Introduc t ion 

Iterative  reconstruction  of  distorted  signals  has  received  much  attention 
in  the  engineering  literature.  Many  algorithms  have  been  presented  for  dif¬ 
ferent  models  of  signals.  The  reader  is  refered  to  [6]  for  a  comprehensive 
review. 

In  this  paper,  we  will  present  an  approach  which  unifies  a  number  of 
important  algorithms  in  the  restoration  of  linearly  distorted  signals.  The 
basic  tool  which  we  will  use  is  that  of  iterative  solution  of  linear  operators 
in  Hilbert  spaces.  The  advantages  of  this  approach,  which  is  based  on  a 
result  given  by  Bialy  [lj ,  are  the  following: 

1.  Several  apparertly  disconnected  algorithms,  some  of  which  have  recently 
received  much  interest,  can  be  considered  special  cases  of  Bialy's  itera¬ 
tion. 

2.  All  these  algorithms  can  be  shown  to  be  convergent  using  a  rather  general 
tool . 

3.  A  simple  generalization  of  the  basic  iterative  procedure  will  be  shown  to 
provide  some  new  restoration  algorithms  which  perform  fast  reconstruction 
of  the  sought  signal. 

Section  II  reviews  some  fundamentals  of  linear  operators  in  Hilbert 
spaces.  Special  emphasis  is  put  on  pseudoirverse  solutions  and  Bialy's  itera¬ 
tion  for  non-negative  symmetric  operators.  In  Section  III,  we  show  that  th  s 
iteration  can  be  used  to  obtain  the  iterative  procedures  mentioned  in  the 
abstract  of  this  paper.  In  particular,  we  obtain  a  generalization  of  the 
Papoul  is-Gerchberg  algorithm  for  the  continuous  extrapolation  problem.  In 
Section  IV,  we  show  how  a  simple  generalization  of  Bialy's  iteration  provides 
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some  useful  recursive  techniques  for  restoration.  Some  numerical  examples 
showing  the  performance  of  these  algorithms  are  presented  in  Section  V,  where 
the  application  problem  is  continuous  band-limited  signal  extrapolation.  A 
numerical  comparison  of  these  algorithms  with  the  Papoulis-Gerchberg  procedure 


is  presented. 


II .  Basic  theory 


Let  as  recall  what  is  meant  by  bounded  and  compact  linear  operators  in 
Hilbert  spaces.  Let  Hj  be  two  Hilbert  spaces  and  A:  H^  ->  H2  a  linear 
operator.  We  say  that  A  is  bounded  (also  continuous)  if  there  exists  a  real 
number  C  such  that 


llA  x||2  ^  c  1 1x1  ^  for  all  x  6  Bj 
where  II  II.  denotes  the  norm  in  Hj. 

The  operator  A  is  called  compact  if  it  maps  every  bounded  set  S  C  onto  a 

set  A(S)  whose  closure  A(S)  is  compact.  In  other  words.  A  is  compact  if  and 
only  if  for  every  bounded  sequence  {x  n  €  N}  C  H.  ,  there  exists  a  subse- 

Tl-ll, 

quence  Un  ,  k  €  N)  and  y  6  H2  such  that  A(xn  )  - >  *  “>  “•  T*®  reader 

k  k 

is  refered  to  [7]  for  further  theoretical  details. 

Obviously,  if  a  linear  operator  is  compact  it  will  also  be  bounded.  The 
converse  does  not  hold  in  general.  However,  if  is  of  finite  dimension  both 
class  of  operators  coincide.  The  adjoint  of  A  is  another  linear  operator  A* : 
**2  ~ ^  Hi  characterized  by  the  following  identity: 


<Aty, 


<y. 


where  <  ,  >g  denotes  the  inner  product.  A  linear  operator  A:  H  ->  H  is  called 
i 

symmetric  if  A  =  A1.  In  that  case,  we  say  that  A  is  non-negative  if  <Ax,x>  2  0 
for  all  x  G  II.  We  concern  ourselves  with  iterative  solutions  to  the  linear 
problem  Ax  =  y,  where  A:  Ci  ->  n2  is  bounded  and  y  6  H2  is  given. 
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Frequently,  it  happens  that  y  does  not  belong  to  the  range  of  A  and 
therefore,  thre  is  no  x:  Ax  =  y.  In  that  case,  one  may  attempt  to  find  the 
minimum  norm  least  squares  solution.  However,  for  infinite-dimension  spaces 
this  approach  is  not  always  successful  because  the  "least  square"  solutions 
may  fail  to  exist.  We  need  to  recall  some  related  results  for  our  applica¬ 
tions. 

It  is  well-known  that  the  range  of  a  bounded  operator  may  not  be  closed. 
The  situation  is  even  worse  for  compact  operators  since  it  can  be  proved  that 
the  range  of  such  an  operator  is  "almost  never"  closed.  Undoubtedly,  this 
result  is  the  main  drawback  for  a  pseudoinverse  approach  to  solving  the  opera¬ 
tor  equation  Ax  =  y,  because  most  of  distortion  equations  in  signal  processing 
are  given  by  compact  operators.  The  following  lemmas,  which  are  proved  in  ref. 
[8],  help  in  understanding  the  matter,  and  will  be  useful  for  the  remainder  of 
our  paper. 

Lemma  1 

Let  A:  ->  H2  be  a  linear  bounded  operator.  For  a  fixed  y  6  H2 »  let  S  = 
(x  6  Hj:  Ax  =  Qy)  and  N  =  (x  €  A*Ax  =  Aly] .  Then  S  =  N.  (Q:  E2  ->  aTh^T 
is  the  projection  operator  onto  the  closure  of  the  range  of  A.)  The  equation 
A*Ax  =  A*y  is  recognized  as  the  normal  equation  for  A.  It  is  obvious  that  if 
Qy  does  not  belong  to  the  range  of  A:  R(A)  then  N  =  S  =  Therefore,  since 
R(A)  may  not  be  closed,  there  exists  many  points  y  €  H2:  Qy  g  R(A) .  In  other 
words,  N  will  not  be  empty  iff  y  €  R(A)  +  R(A)'*'  (1  denotes  orthogonal  set). 
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T  1emng  2 

For  a  fixed  y  6  H^,  the  set  of  least  squares  solutions 

{u  6  H1 ;  J |Au  -  yll  *  inf  tilAx  -  yll,  x  6  H^} } 

coincides  with  the  set  of  solutions  of  the  normal  equation  A*Au  *  Afcy. 

From  the  comments  given  above,  it  is  clear  that  the  set  of  least  squares 
solutions  will  not  be  empty  if  and  only  if  y  6  R(A)  +  R(A)1.  In  that  case, 
this  set  will  be  closed  and  convex  and  therefore,  there  will  be  an  element  u  - 
A+y  which  has  minimum  norm  among  all  which  satisfy  A*Au  =  Aty. 

Another  simple  but  very  important  property  is  the  following: 

Lemma  3. 

If  y  €  R(A) ,  then  A+y  is  the  minimum  norm  solution  of  the  linear  equation 
Ax  =  y. 

Lemma  3  says  that  if  a  solution  to  the  problem  Ax  =  y  exists,  then  the 
minimum  norm  solution  will  make  sense  and  will  coincide  with  the  generalized 
inverse  A+y.  Ibis  is  a  very  simple  consequence  of  the  fact  that  y  6  R(A) 
ensures  that  the  normal  equation  A*Ax  =  A*y  has  the  same  set  of  solutions  of 
Ax  =  y. 

One  would  like  to  have  pseudoinverse  solutions  for  every  y  €  .  How¬ 

ever,  as  we  have  shown  above  this  will  be  possible  iff  R(A)  is  closed.  In  that 
case,  the  generalized  inverse  A+:  H2  ->  Hx  is  a  well-defined  bounded  operator. 
The  boundeness  of  A+  shows  that  finding  a  pseudosolution  A+y  is  a  stable  prob¬ 
lem,  i.e.  small  perturbations  in  the  data  will  produce  small  changes  in  the 
pseudosolution  A+y.  As  we  have  mentioned  above,  this  will  be  "almost  never" 
the  case  if  A  is  compact: 


Lemma  4. 

If  A:  ->  H2  is  compact  and  R(A)  is  closed  then  A  is  degenerate,  i.e. 

R(A)  is  of  finite  dimension. 

Some  examples  of  compact  operators  may  clarify  the  matter.  Let  ns  suppose 
that  oar  distortion  can  be  written  either  as 

A  a 

g(y)  *  /  h(x,y) f (x)dx,  y  6  (-b,b)  (for  all  f:  7  |f(x)l2dx  <  ®) 

-a  -a 

for  a  continuous  model  or  as 

g(n)  =  £  h(m,n)f(n),  m  €  Z  (for  all  f:  £  lf(n)|2  <  <*») 

n  €  Z  n  €  Z 

for  a  discrete  model. 

In  both  cases,  under  rather  general  conditions  on  h  it  can  be  shown  that 
the  corresponding  distortion  operator  is  compact.  Sometimes,  the  situation  is 
even  worse  because  the  range  is  not  only  non-closed  but  also  dense  (i.e.  R(A) 
-  Hj)  .  Practical  terms,  this  means  that  if  the  given  data  g  is  contaminated 
with  some  additive  noise  ij,  the  problem  becomes  intractable  from  a  generalized 
inverse  point  of  view.  This  is  because  R(A)^  =  {0}  and  therefore,  A+g  will 
never  exist  if  the  noise  has  any  component  which  is  outside  of  R(A)  (this  is 
almost  always  the  case).  An  example  of  dense  range  is  provided  by  the  set  of 
band-limited  functions,  with  a  fixed  bandwidth  D.  It  is  well  known  that  this 
set  is  dense  in  the  set  of  finite  energy  functions  over  an  interval  [-a, a] 
(see  [9], [2]).  We  hope  to  shed  more  light  on  this  problem  in  section  III. 2. 

In  what  follows  we  will  state  the  Bialy  iteration  which  is  also  useful 
for  computing  generalized  inverse  solutions  of  Ax  =  y.  This  iteration  is  the 
main  core  of  the  next  section,  and  provides  the  basic  tool  for  the  announced 


unification  of  algorithms. 

To  this  end,  if  A  is  a  bounded  linear  operator,  ve  denote  by  | | A | |  the 
inf imum  of  the  numbers  c:  llAxll  i  clixll,  for  all  x  6  .  Ve  also  denote  by 

P  the  orthogonal  projection  onto  the  kernel  of  A:  Eer(A)  *  (x  6  H,  :  Ax  ■  0} . 

Theorem  1  (Bialy,  [1]) 

Let  A:  H  ->  H  be  a  linear  bounded  non-negative  operator.  For  y  e  H,  xq  € 
H  consider  the  iterative  process 

xn+l  "  *n  +  °<y  “  ^n-l*  <*> 

where  0  <  a  <  2/ I  I A I  I . 

Then,  the  sequence  (xn>  n  2.  0]  converges  if  and  only  if  Ax  =  y  has  a 

solution.  In  that  case  xq  ->  px<j  +  i,  where  x  is  the  minimum  norm  solution. 

n->« 

Ve  would  like  to  make  some  remarks  about  Theorem  1.  It  is  clear  that  if 
the  initial  approximation  xq  is  zero  then  {xnI  will  approach  the  minimum  norm 
solution  of  the  equation  Ax  -  y.  The  theorem  also  says  that  this  will  happen 
iff  the  equation  has  at  least  one  solution. 

Judging  from  the  appearance  of  (1)  it  may  be  said  that  recursion  (1) 
tries  to  compute  a  fixed  point  of  the  mapping  Gx  *  ay  +  (I  -  aA)x.  However,  it 
cannot  be  said  that  the  fixed  point  and/or  the  iterative  procedure  make  sense 
because  of  a  contractive  property  of  G.  In  fact,  this  situation  will  almost 
never  occur.  The  reason  of  this  assertion  is  given  by  the  following. 

Lemma  5 

If  A:  H  ->  H  is  a  bounded  linear  operator  such  that  iteration  (1)  con- 
V€r8es  *11  y  €  H  for  some  x  and  a  0  then  A  cannot  be  compact,  unless  H 
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is  of  finite  dimension. 

The  proof  of  this  lemma  is  relegated  to  the  reader.  Lemma  5  says  that  if 
A  is  a  compact  operator  and  the  dimension  of  H  is  not  finite  (and  therefore, 
including  most  of  the  cases  we  are  interested  in)  recursion  (1)  must  be  diver¬ 
gent  for  some  y.  In  particular,  I  -  XA  will  not  be  a  contraction  mapping 
irrespective  of  the  choice  of  A. 

A  relevant  characteristic  of  the  hypotheses  of  Theorem  1  is  that  A  is 
assumed  to  be  non-negative,  excluding  apparently  many  operators  for  which  this 
condition  is  not  met.  However,  Bialy's  theorem  can  be  used  to  compute  itera¬ 
tively  the  minimum  norm  least  squares  solution  of  any  bounded  linear  operator. 
This  result  will  be  obtained  very  easily  if  we  recall  that  the  minimum  norm 
least  squares  solution  (whenever  exists)  is  the  minimum  norm  solution  of  the 
normal  equation  A*Ax  *  A*y.  Then,  Bialy's  theorem  can  be  applied  because  A*A 
is  a  non-negative  linear  bounded  operator.  Thus,  we  have: 

Theorem  2 

Let  A:  ->  H2  be  a  bounded  linear  operator,  y  6  R(A)  +  R(A)"^,  consider 

the  iterative  equation: 

Xo  "  0 

*n  =  xn-l  +  aAt(y  "  Axn-1>'  n  2  1  (2) 

whe’-e  0  <  a  <  2/||AtA||.  Then  t*n)  converges  to  the  minimum  norm  least  s- uares 
solution  A+y. 

It  is  worth  noticing  that  Theorem  2  assumes  y  €  R(A)  +  R(A)'* 
fore,  the  sought  generalized  inverse  solution  A+y  exists. 


and  there- 
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To  end  this  section,  we  would  like  to  point  out  some  results  which  ere 
connected  to  those  of  Theorem  1  end  Theorem  2.  Two  special  esses  of  Bialy's 
theorem  were  proved  esrlier  for  integrel  operators,  which  are  e  case  of  com¬ 
pact  operators  arising  very  often  in  practical  applications.  In  1956,  Fridman 

a 

([10])  proved  Theorem  1  for  the  case:  A  is  given  by  Af(x)  *  /_  h(x, t) f (t)dt , 
f  is  any  finite  energy  function,  i.e.  f  €  L  (-a,a),  and  the  kernel  h  is  posi¬ 
tive  (and  symmetric:  h(x,t)  *  h(t,x)) .  In  1951,  Landweber  ([11])  proved 

a 

Theorem  2  for  the  case:  Af(x)  =  h(x,t)f(t),  removing  the  assumptions  made 
on  h(x,t).  In  both  cases,  h(x,t)  must  define  a  compact  operator,  a  condition 
that  is  often  met. 

A  final  remark  is  in  order;  it  can  be  easily  obtained  from  Theorem  2  and 
the  discussion  on  pseudosolutions  presented  in  this  section,  that  iteration 
(2)  will  appraoch  the  minimum  norm  solution  of  the  equation  Ax  =  y,  whenever  y 
€  R(A) .  To  see  that  we  just  recall  that  [x  6  H^:  AtAx  =  Aly]  =  {x  €  E^;  Ax  = 
y]  for  y  €  R(A) .  In  particular,  if  the  solution  exists  and  is  unique  it  will 
be  also  obtained  by  the  procedure  (2). 
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III.  Applications 

In  this  section  we  will  show  several  applications  of  the  resnlts  dis¬ 
cussed  in  Section  II. 

III.l  Van  Cittert-type  algorithms. 

2 

Ve  now  consider  a  continuons-continnous  deconvolution  problem.  Let  L  (B) 
denote  the  Hilbert  space  of  finite  energy  functions  defined  on  B,  i.e.  L  (B) 
=  {f:  B  ->  R:  /g  |f(t)|2dt  <  *}  .  Let  h  be  a  function  such  that  the  following 
linear  operator  is  bounded: 

A:  L2(S)  ->  L2(T) 

f  ->  Af ( t )  =  /g  h ( t-s) f ( s)ds ,  t  G  T 

(If  T  or  S  is  bounded,  and  h  satisfies  Aj/jJh (s-t)  |2dsdt  <  ®  then  A  will  be 
bounded.  In  that  case,  it  can  be  proved  that 

,|AfN  2  1  {/s/T|h(s-t)|2dsdt)1/2  I  If  I  I 

L2(T)  a  A  L2(S) 

where  1 1 y  1 1  stands  for  the  norm  (/n ly ( t ) 1 2dt } 2 .  Another  case  for  which 

L2(B) 

A  is  bounded  will  be  obtained  if  the  function  h  has  compact  support,  that  is 
to  say,  h(s)  =  0  if  s  6  C  where  C  is  a  compact  set  in  Rn.) 

If  S  =  T,  J,,  denotes  truncation  to  S  and  h  satisfies  the  additional  pro¬ 
perties  /  |h(t)|2dt  <  ®  and  h(<a)  2  0  for  all  w  ®  Rn,  where  h  denotes  Fourier 

Rn 

transform,  then  A  is  a  non-negative  operator.  To  see  this, 

<Af,f>  =  /o(Af) (s) .f (s)  =  /  (h*Jsf) (s)Jsf (s)ds 

L2 ( S)  b  Rn  b  b 

where  the  symbol  *  stands  for  convolution,  defined  over  Rn.  By  means  of 

Parseva  1  ’  s  equality,  we  obtain  <Af,f>  =  /  (Hwj^fJ^w)  {  ( J  „  f  )A]  (w)dw.  Rut 

L2(S)  Rn  b  6 


/  h(w)  .1  (Jsf)A(w)  |2dw,  which 

Rn  S 


(k*Jcf)(w)  *  h(w) . (J«f  )A  (w) ,  then  <Af,f>  - 
S  S  L2 (S) 

is  always  non-negative. 

We  now  assume  that  g(s),  s  6  S  is  the  output  of  our  continuous  system 
defined  by  A.  If  we  are  interested  in  recovering  the  input  f(s).  s  €  S,  we  can 
apply  Theorem  1  to  obtain  a  sequence  {fQ}  given  by 

£o  *  0 

£n+l  ‘  f»  +  “<*  -  <3> 

which  converges  to  the  minimum  energy  signal  that  produces  the  output  g. 
Another  way  of  writing  equation  (3)  is 

^n+i ( *)  *  +  ~  /gh(s-t) fn(t)dt) ,  s  €  S, 

or  equivalently,  fn+l  (s)  *  ag(s)  +  (fn(s)  "  o/sh(s-t)fn(t)dt) . 

This  is  a  Van  Cittert-type  recursion  whose  convergence  is  ensured. 
Several  remarks  are  in  order.  Perhaps  the  most  important  observation  is  that  h 
may  have  zero  frequencies  without  affecting  the  convergence  of  the  procedure. 
It  is  also  clear  that  many  choices  of  a  can  be  tried  whenever  0  <  a  <  2/ | I A I  I 
(if  S  is  bounded  we  can  chose  any  a  which  satisfies 

0  <  a  <  2/(/s/s!h(*-t) |2dxdt)1/2)  . 

The  classical  Van  Cittert's  algorithm  is  for  the  case  S  =  T  =  Rn.  It  is 

this  assumption  what  makes  the  proof  of  Van  Cittert's  iteration  (3)  so  simple 

if  a  is  chosen  to  satisfy  |l  -  ah(w)l  <  1  whenever  h(w)  #  0  ([6]).  Therefore, 

A 

i  f  S  R  ,  under  the  more  stringent  condition  h(w)  >_  0,  Bialy's  iterati-'  pro¬ 
vides  a  non-trivial  extension  of  the  classical  version  of  Van  Cittert's  algo¬ 
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We  next  consider  a  discrete-discrete  deconvolution  problem.  In  this  case, 
the  underlying  Hilbert  space  is  1-(B)  -  {a  :  m  €  B,  7  la  1^  <  *)  »  where  B 

m  e  b 

is  a  subset  of  Zn.  Let  h(m)  be  a  sequence  such  that  the  operator 

A:  12(B)  ->  12(C) 

x(m),  m  6  Zn  ->  £  h(k-m)x(m),  k  6  C 

m  e  B 

is  bounded.  Several  conditions  on  h,  similar  to  those  given  for  the  continuous 
case,  can  be  found  to  ensure  the  boundness  of  A. 

If  h  satisfies  £  |h{m)l^  <  eo,  B  =  C  and  the  Fourier  series  of 

m  €  Zn 

h:  7  h (m) e^rrmu)  1  0  for  all  <*>»  then  Bialy’s  theorem  will  apply.  Thus,  the 

m  €  Zn 
iteration 

*£k)  =  x^j  +  afgttJ-b'JgX^tk))  ,  k  6  B  (4) 

will  converge  to  the  minimum  norm  solution  of  the  problem:  g (m)  =  (h*Jgf) (m)  , 
m  €  B,  provided  that  at  least  one  solution  exists.  Equation  (4)  can  be  written 
as  follows 


xm(k)  ag(k)  +  xm_1(k)  _  Q  y  h(k_j )Xm_1 ( j ) ,  k  €  B 

j  E  B 


(4') 


with  0  <  a  <  2/ I  I A 1  I .  A  simple  rationale  for  choosing  a  is  0  <  a  <  2/ ( sup  I b I ) 


Equation  (4)  (or  its  equivalent  form  (4’))  is  a  Van-C i t t e r t ’ s  recursive 
formula,  when  the  model  for  the  observed  and  unknown  signals  are  both 
discrete. 


It  is  worth  noticing  that  equation  (4')  can  also  be  written  by  using  an 


operator-type  notation: 
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x^k)  =  ag  (k)  +  {(6  -  ohJ^Jgx^}  (k)  .  k  6  B  ,  (5) 

where  6:  Zn  ->  C:  6(0)  *  1.  6(k)  =0  if  k  M. 

Recursion  (5)  was  also  considered  in  [6]  for  the  special  case  in  which  B 

is  bounded.  Under  this  assumption,  equation  (S)  was  shown  to  converge  in  ref. 

[12].  However,  in  [12]  it  was  also  proved  that  if  B  is  bounded  a  can  be  chosen 

independently  of  B  so  that  (6  -  ah)*Jg  is  a  contraction  mapping  under  mild 

A 

assumptions  on  h  (which  include  the  case  h(ui)  2  0  for  all  w  6  R  ).  We  think 
that  it  is  now  understood  what  is  the  theoretical  importance  of  the  a  priori 
constraint  that  the  sought  signal  is  limited  to  the  set  B.  From  lemma  S  and 
related  discussions,  it  is  seen  that  (6  -  Ah)*Jg  will  be  a  contraction  mapping 
for  a  certain  X  and  for  a  rather  general  h  only  if  B  is  bounded.  On  the  other 
hand,  if  B  is  not  bounded,  equation  (S'  was  shown  to  converge  to  the  minimum 
norm  solution  of  the  deconvolution  problem  (Theorem  1),  but  the  contractive 
property  of  (6  -  Xh)*Jg  will  not  hold  in  general. 

To  conclude  this  subsection  we  would  like  to  emphasize  that  for  the 
continuous-continuous  model,  if  the  set  S  where  the  input  signal  f  is  not  zero 
is  bounded,  (I  -  ah)*J<,  will  not  be  a  contraction  mapping  in  general.  This  is 
a  major  difference  between  the  continuous  and  discrete  model,  ^(B)  is  of  fin- 

2 

ite  dimension  if  B  is  bounded  whereas  L  (S)  does  not  have  this  property. 

III. 2  Pse udo inve rse  regularization 

The  deconvolution  problem  that  was  discussed  in  Section  III.l  usually 
requires  a  more  involved  solution  due  to  the  following  facts: 

1.  g  is  given  with  noise  and  therefore  the  solution  to  the  problem  g  =  Af 


may  not  exist. 
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2.  h  may  not  satisfy  h(w)  2  0  lor  all  w. 

3.  The  period  of  time  where  the  observation  g  is  given:  C  may  not  coincide 
with  the  support  6  of  the  sought  signal. 

A  full  answer  to  problems  2  and  3  and  a  partial  solution  to  1  will  be  given  in 
this  section.  To  this  end,  we  will  show  the  convergence  of  an  iterative  recon¬ 
struction  algorithm.  We  will  consider  the  discrete-discrete  model  only.  Simi¬ 
lar  comments  and  results  hold  for  the  continuous-continuous  case.  With  the 
same  notation  as  in  Section  III.l,  let  us  suppose  that  our  observation  g(m),  m 
€  C  is  given  with  noise.  Assume  that  g  €  R(A)  +  R(A)^.  Then,  the  following 
problem  will  always  have  a  solution:  AfcAf  =  Alg  where  A*  denotes  the  adjoint 
of  A  (see  Section  II).  For  the  convolution  case: 

(Af) (m)  =  7  h(m-k) f (k) , 

k  6  B 

m  €  C  then,  (Alq) (k)  =  £  h (-k+m) q (®) >  k  £  B.  This  means  that  A1  is  also 

m  e  c  _ 

given  in  terms  of  a  convolution  where  the  new  kernel  is  h  (-m) .  Specifically, 
ArA  is  given  by 

(AtAf)(j)  = 7(7  h (m-j ) h (m-k) ) f (k) ,  j  €B 
k  e  B  m  €  C 

which  is  always  non-negative,  as  it  was  pointed  out  in  Section  2. 

Thus,  we  can  apply  Theorem  2  for  computing  the  minimum  norm  least  squares 
solution  of  Af  =  g  ,  if  g  G  R(A)  +  R(A)*,  by  means  of  the  iteration: 

Xo  =  0 

Xm  =  xm-l  +  «At(g  “  Axm-1>  "  >  1  (6) 

where  0  <  a  '  2/ ( |  I  A  I  J  2 )  . 
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An  equivalent  expression  for  (6)  is  obtained  by  replacing  A  and  A 1 

*0  =  0 

k  6  B,  x  (k)  =  a  E  h(m-k)g(m)  +  x  1(k)  -  a  E  (  E  h(m-k)h(m-j) ) f ( j) 
m  €  C  j  6  B  m  €  C 

(o'') 

It  is  worth  pointing  out  that  equation  (6)  or  (6')  will  not  converge  if  the 
noisy  data  g  6  R(A)  +  R(A)^.  However,  we  think  that  this  approach  is  useful 
for  understanding  the  limitations  of  the  technique  and  for  setting  a  condition 
to  ensure  convergence  or  divergence  of  the  iteration. 

A  particular  case  is  obtained  for  B  =  C  =  Zn.  In  that  case,  the  technique 
that  consists  of  convolving  the  equation  h  •  f  =  g  with  h(-m)  has  been  pro¬ 
posed  independently  by  several  authors  ([6], [13])  but  the  approaches  used  were 
conceptually  different.  For  B  =  C  =  Zn,  the  operator  A*A  is  given  by  a  convo¬ 
lution  whose  kernel  is  h(m)  *  h(-m).  Then,  the  transfer  function  of  the  system 
AtA  is  lh(w)|^,  w  6  Rn.  Since  ih(w)|^  is  always  non-negative,  Van-Cit tert '  s 
algorithm  applies  to  the  equation  [h ( -m) *h (m) ] *f  =  h(-m)*g,  obtaining  the  fol¬ 
lowing  frequency-space  recursion: 

*Q(w)  =  0 

*m(w)  =  afi(w)  g  (w)  +  (1  -  a|  h  (w)  I  2)  xjd_1  (w)  (7) 

is,  of  course,  equivalent  to  (6’)  wh.  B  =  C  =  Zn.  It  is  worth  pointing 
j  *.  there  should  be  a  solution  to  the  problem  h(m)*f  =  g  in  order  to 
,  :  :  n i t tf-ene rgy  discrete  signal  whose  Fourier  transform  is  the  limit  of 


ir»)ve,  the  pseudo  inverse  approach  provides  a  full  answer  to 
jnd  C  are  any  subset  of  Zn ,  and  the  convergence  of  (6)  is 
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characterized  in  terms  of  the  data  g. 

To  end  this  subsection  we  would  like  to  remark  that  if  g  0  R(A)  +  R(A)^ 

then  (6)  will  not  converge.  Therefore,  caution  is  recommended  when  the  itera¬ 
tion  (6)  (or  (6'))  is  used  for  theoretical  derivations.  On  the  other  hand, 
when  (6)  (or  (6'))  is  implemented  numerically,  a  finite  piece  of  the  signals 
is  used  and  therefore  convergence  of  the  iteration  is  guaranteed.  The  concep¬ 
tual  point  is  that  for  implementation  purposes,  the  underlying  model  for  the 
distortion  is  g(m)  =  (h*Jgf) (m) ,  m  £  C  where  both  B  and  C  are  finite.  Then, 
the  pseudoinverse  solution  will  always  exist  and  can  be  approximated  by  means 
of  (6  ’)  . 

III. 3  Papoul i s-Gerchberg  '  s  iteration 

Let  us  assume  that  g:  F  ->  C  is  a  piece  of  a  O-band-limited  function 

(i.e.  g(w)  =  0,  w  6  0)  where  F  M  is  an  open  subset  of  Rn.  Let  us  suppose 

that  the  complete  function  g  satisfies  the  finite-energy  constraint: 

/  Ig  (x)  |2dx  <  <». 

Rn 

Papoulis  [2]  and  Gerchberg  [3]  proposed  the  following  algorithm  for  com¬ 
puting  the  continuation  of  g  or  its  Fourier  transform: 

®o  =  0 

gm  =  sincQ*(JFg  +  (I-JF)gm_1) ,  m  2  1  (8) 

where  sinc^  denotes  the  function  whose  Fourier  transform  is  the  indicator  of 

Q. 

In  ref.  [2],  equation  (8)  was  shown  to  be  convergent  to  g  in  the  energy 
norm  for  the  one-dimensional  case.  In  ref.  [14],  another  approach  was  shown  to 
prove  convergence  of  (8)  which  is  also  valid  for  the  mul t i -di men s i ona 1  case. 
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However,  in  [15],  this  algorithm  was  presented  as  a  special  case  of 
Landweber's  iteration  ([11]).  The  underlying  operator  equation  is  (Af) (x)  * 

g(x)  where  A  is  an  integral  operator  given  by 

(Af) (x)  =  /  f(w)e_2niwxdw,  x  6  F  (9) 

Q 

f  €  L  (G)  and  G  is  assumed  to  be  bounded.  It  is  obvious  that  the  sought  solu¬ 
tion  is  f  =  g  and  is  also  unique.  We  can  now  apply  Theorem  2  to  get  a  recur- 


fo  =  0 

fm  =  fm-l  +  oAt(*  ‘  A£m-1>  .  <10) 

where  0  <  a  <  2/(||AtA||)  and  A*  is  the  adjoint  integral  operator  given  by 
(Ath)(w)  =  /p h(x)e2niwxdx.  w  €  Q. 

It  is  very  easy  to  verify  that  I  I A A I  I  <,  1;  then,  a  =  1  i  s  an  admissible 
value  and  from  (10) 


fo  "  0 


fm  *  Vl  +  A*(g  -  Af.^) 


(10') 


will  converge  in  the  energy  norm  to  the  unique  solution  f  of  the  equation 

v  v 

(Af)(x)  =  g(x),  x  G  F.  But  in  that  case,  f  ->  f  also  in  the  energy  norm  (v 

m_>”  v 

denoi  inverse  Fourier  transform).  Since  f  f  are  supported  on  C,  fn(x)  =  / 

v  n  G 

i  ( ;,)  p.-"  T*xwdw  and  f(x)  =  /  f  (w)  e_2;Tiwxdw .  x  6  Rn.  We  now  apply  inverse 

'*  G 

<ansform  to  both  sides  of  (10’): 


fo  =  0 


fm  =  fm-i  +  («(x) 


/  fm_,(w)e"2lTixWdw)e 


o  -  ,  V  ,  V 

i',1I-dx} 
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If  we  cell  g  »  f  in  >  0  we  will  obtain  the  recursion 

IQ  ID  — 


*o  -  0 


gm  "  8m-l  +  sincQ*<JF«  -  Jpg^) 


(11 ') 


Since  gm_j  is  Q-band-limited,  equation  (11')  is  equivalent  to  the  following: 


g  =  0 


:Q*(JFg  +  (I-JpJg^).  m  >  1 


Bm  ■  sincfl* 


Equation  (12)  is  the  Papoulis-Gerchberg  iteration  (8) . 

We  will  now  show  a  generalization  of  (8)  to  cases  where  the  low-pass  con¬ 
volution  is  performed  by  some  other  operator. 

To  this  end,  we  need  to  assume  some  further  information  related  to  g  =  f. 

A 

Let  us  suppose  that  for  certain  non-negative  bounded  function  h(w) ,  w  €  D, 

A 

h(w)  =  0,  w  (B  Q,  g  satisfies: 


Ig  (w)  |- 


/  -  dw  <  - 

h(w) 


2  2 

Then,  if  we  consider  the  operator  A:  L  (0)  ->  L  (F) : 

(As)  (x)  =  fQ  h1/2U)  .e2jliwxs(w)dw,  x  E  F  (14) 

2 

the  equation  g(x)  =  (As)(x),  x  E  F  will  have  a  solution  in  L  (Q)  (which  is 
obviously  g/h1^2) .  It  readily  follows  that  the  solution  is  also  unique.  We  can 
now  apply  Theorem  2  to  the  equation  g(x)  =  (As) (x) ,  x  E  F  for  A  given  in  (14) 
to  get  an  iterative  procedure: 


fo  =  0 


fm  =  fm-l  +  «At(8  ~  .  m  >  1 
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which  converges  to  the  unique  solution.  More  specifically  equation  (15) 
becomes 

fo-° 

f  (w)  =  fm_i(w)  +  o/fi1/2(w)e2niwx(g(x)  -  /£1/,2(z)e“2nixxf  .  (z)dz)dx 

F  D 

(15') 

If  we  now  multipy  both  sides  of  (15')  by  £^2(w)e“2n*wy,  integrate  respect  to 
w  €  Q,  and  call  g_(y)  =  /  f  (w)£1/,2(») e-2iriwydw,  we  will  obtain 

m  q  ® 

«o  =  0 

y  €  Rn,  gm(y)  =  gm_1(y)  +  a  /  {/  ft(w)e2niw^x  y^dw](g(x)  -  gffl_1(x))dx  (16) 

F  Q 

If  we  call  h(z)  =  /  h (w) e2niwzdw,  equation  (16)  will  become 
g0  =  0 

*m  =  *m-l  +  *h(-z)*JF(S  -  gm_1)l  m  >  1  (17) 

which  converges  to  g  uniformly  over  compact  sets  in  Rn,  when  0  <  a  < 

2/ (sup |h (w) | )  and  (13)  is  satisfied. 

w 

Two  well-known  discrete  algorithms  for  extrapolation  can  be  considered  a 
sampled  version  of  equations  (8)  and  (11')  (see  [16], [17]).  In  addition,  some 
new  algorithms  for  solving  the  discrete  extrapolation  problem  ([18])  can  also 
be  interpreted  as  a  sampled  approximation  of  equation  (17). 

III. 4  Iterative  extrapolation  of  ’nf  inite-extent  discrete  signals. 

Let  F  oe  a  finite  subset  of  Zn  and  z(m),  m  6  F  a  sequence  of  numbers. 
The  discrete  band-limited  extrapolation  problem  consists  of  finding  an  infin¬ 
ite  sequence  y(m) ,  m  6  Zn  such  that  y(m)  =  z(m),  m  €  F  and  y(m)  is  O-band- 
limited,  i.e.  y(w)  =  y  y(m)e-2,'I*n:w  =  0  if  w  6  Q  (a  fixed  bounded  set  of 

m  c  7^1 


frequencies  QC  [— 1 .1] n  (see  [4] , [16] , [17] , [18] ) . 


The  solution  to  this  problem  is  non-unique  ([6], 14]).  In  ref.  [4]  it  was 
shown  that  the  minimum-norm  discrete  extrapolation  y  can  be  computed  by  means 
of  the  following  two-step  procedure: 

1.  solve  for  x:  I  sinc0(k-m)x(m)  =  z(k),  k  €  F  (18a) 

m  €  F 

2.  compute  £  sine (k-m) x (m)  =  y(k),  k  6  Zn  (18b) 

m  €  F 

Then,  it  was  shown  that  y  can  be  computed  by  the  following  iterative  algo¬ 
rithm: 


y0  -  o 

^  €  Zn*  =  ym-l^  +  °  ^  sincc(k-j)(z(j)-y  .(j)),  m  2  1  (19) 

j  €  F 

for  0  <  a  <  2.  (Both  results  were  extended  for  arbitrary  multidimensional  F 
and  G  in  [18];  for  a  relationship  between  this  discrete  solution  and  the  con¬ 
tinuous  extrapolation  problem  given  in  III. 3,  see  [18].) 

Perhaps,  the  earliest  reference  to  the  technique  given  by  (18a)-(18b)  is 
Yao  [19]  who  addressed  this  problem  under  a  rather  different  name  and  by  using 
a  quite  general  approach. 

The  fact  that  iteration  (19)  computes  the  same  sequence  as  that  of 
(18a)-(18b)  is  very  simple.  In  this  section,  we  will  show  that  the  iteration 
(19)  can  be  obtained  from  Bialy's  iteration  for  a  certain  operator  equation 
problem.  The  minimum  norm  property  of  the  limit  sequence  will  be  readily 
derived  as  a  byproduct. 

Let  A:  L  (G)  ->  ^(F)  be  the  following  linear  operator: 


(Af)(m)  -  /  f(w)e2ltinwdw,  m  6  F 
0 

9 

It  is  clear  that  A  is  bounded  when  L  (0)  and  l^F)  are  equipped  with  the  norms 

•f|f(w)|2dw  and  7  |z(m)|2  respectively.  It  is  clear  that  the  discrete  extra- 

12  m  €  F 

polation  problem  can  be  put  in  this  equivalent  way: 

find  f  G  ( £2) :  (Af)(m)  =  z(m),  m  €  F  (20) 

From  Parseval's  formula,  it  is  seen  that  the  minimum  norm  extrapolation 
corresponds  to  minimizing  llfll^  where  f  satisfies  (20).  We  can  now  solve  (20) 
by  means  of  Bialy's  iteration.  To  this  end,  we  need  to  compute  A*.  It  is  sim¬ 
ple  to  verify  that,  if  s  €  ^(F)  then 

(Afcs) (w)  =  7  s<m)e"2nim\  w  G  [-l.l]n. 

m  €  F 

Thus,  Bialy's  iteration  given  by  Theorem  2  becomes 

f0  =  0  (21) 

f_(w)  =  f  ,(w)  +  a  Z  (z(k)  -  /  f_  ,  (z)  e2nizkdz )  e-2nikw,  w  G  0. 

k  G  F  0 

9 

and  f  converges  to  the  minimum  norm  solution  of  (20)  in  the  L  (0)  norm. 

ZD 

Therefore,  /„ f  (w)  e27I*wkdw,k  G  Zn  approaches  the  minimum  norm  0-band-limited 
extrapolation  y(k),  k  G  Zn  when  m  ->  ®  in  the  ^(Z0)  norm  (Parseval's  for¬ 
mula).  Then,  if  we  call  y  (k)  =  /nf  ( w)  e2rt*kwdw ,  k  G  Zn,  equation  (21)  becomes 

ID  Ai  ZD 

yo  =  0 

k  6  Z  ’  ym^  =  ym-l(k'  +  Q  .  £  sincQ(k-j)  ( z  ( j  )  ( j  )  )  ,  m  >.  1  (19) 

j  G  F 


and  convergence  to  y(k) ,  k  G  Zn  is  ensured  for  0  <  a  <  2. 
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A  final  remark  is  in  order.  The  operator  A  given  by  (20)  satisfies  R(A)  = 
^(F)  because  F  is  finite,  therefore  iteration  (21)  or  equivalently  (19)  will 
always  converge  to  the  minimum-energy  solution  of  the  problem.  This  means  that 
the  algorithm  does  not  distinguish  signal  from  noise. 

III. 5  Iterative  extrapolation  of  periodic  discrete  signals. 

Another  related  discrete  approach  to  band-limited  extrapolation  is  to 
solve  the  following  problem: 


Given  z(k),  k  =  -k0,...,k0  <  N 
Find  y (k) ,  -N  <  k  <  N: 


y(k)  =  z(k) .  -k<  k  <  k 
o  —  —  o 


I  y(k)e-2nikn/M  =  0.  In!  >  k 
k=-N  ( 


where  M  =  2N  +  1 . 


In  this  case,  the  band-limited  property  of  y(k) ,  -N  (  k  (  N  is  given  in  terms 

of  the  discrete  Fourier  transform  (DFT) :  £  y (k) e~2nikn/M . 

k=-N 

In  ref.  [5]  the  following  iterative  algorithm  for  computing  the  extrapo¬ 
lation  (22)  was  shown  to  be  convergent: 


y0(k)  =  0,  -N  i  k  (  N 


rPn(k).  -k0  i  k  i  kc 


yn  =  I DFT  i 


C,  otherwise 


(23a) 


where  B  =  DFT 
n 


r z (k ) ,  -k  <  k  <  k 
o  -  -  o 


yn(k) .  ikl  >  ko 


(23b) 


(IDFT  stands  for  the  inverse  discrete  Foarier  transform  given  by 

1/M  *  7  x(k)e2*ikn/M.  n  -  -N . N)  . 

k*-N 

It  is  clear  that  procedure  (23)  incorporates  at  every  iteration  the 
information  available  in  both  time  and  frequency  domains.  In  ref.  [5].  the 
proof  of  the  convergence  of  this  recursion  vas  done  by  means  of  a  certain 
nonexpansive  property  of  an  operator  in  C^. 

In  this  section,  we  show  that  (23)  can  be  also  considered  a  special  case 
of  Bialy's  theorem.  Perhaps,  this  is  the  simplest  of  the  examples  presented  in 
this  section  because  of  the  finite  dimensional  nature  of  the  underlying  Hil¬ 
bert  spaces. 

2kc+l  2k0+l 

Specifically,  let  A:  C  ->  C  given  by  the  IDFT  operator: 

(Ax) (n)  *  £  7°  x(k)e2nikn/M»  -k  £  n  £  k„ 

h  k=-k  °  ° 

o 

It  is  obvious  that  problem  (22)  can  be  restated  as  that  of  finding  a  vector  x 
2k0+1 

€  C  such  that  (Ax) (k)  =  z(k),  -kQ  £  k  £  kQ.  It  is  known  that  this  system 

of  equations  has  a  unique  solution.  We  can  apply  Bialy's  iteration  (2)  for 
computing  the  solution  x.  So,  we  obtain 

xQ(k)  *  0.  -kQ  <  k  i  k0 

Xn  =  xn-l  +  a  At(z  "  ^n-l*'  n  ^  1  (24) 

where  a  can  be  chosen  as  1.  (Here,  A1  is  transpose-conjugate  of  A.) 

We  now  take  M-length  IDFT  on  both  sides  of  (24)  to  obtain 
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yn(k) 


ko 

|rn-l(k>  +  M  h=£°  e2nihk/M  £o  e-2nihm/M  (z  (»)-y  (m) )  (25) 

o  m=-k 

-N  <  k  i  N, 


It  is  easy  to  verify  that  (25)  can  also  be  written  in  the  following  way; 


yo  =  0 


k  N 

y  (k)  =  -  T°  e2nihk/M  r  e-2nihffl/M  (T  ,  .  /T  T  *  \\i  \ 

Jn 'K'  M  e  I  e  z  +  'I-Jk  )yn_i))(m) 

k— — k  m=-N  o  o 


(26) 


where  denotes  truncation  to  £  — kQ , kQ ]  . 

o 

It  turns  out  that  recursion  (26)  is  the  same  as  (23a)-(23b)  and  there¬ 
fore,  the  convergence  of  y^  to  the  sought  extrapolation  is  ensured. 

In  the  derivation  presented  above  it  was  assumed  to  simplify  notation 
that  the  length  of  th?  DFT  is  odd;  2N  +  1 . 


The  advantage  of  this  approach  to  interpreting  iteration  (23)  is  that  it 
is  possible  to  characterize  the  convergence  of  a  similar  procedure  when  the 
number  of  samples  in  the  time  and  frequency  domains  is  not  the  same.  In  such  a 
case,  it  is  obvious  that  the  extrapolation  problem  has  no  solutions  or  an 
infinite  number  of  solutions.  In  both  cases,  the  corresponding  equation  (26) 


will  provide  the  minimum  norm  least  squares  extrapolation. 
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Since  A*A  is  another  linear  non-negative  compact  operator  there  exists  a  fam¬ 
ily  of  eigenvalue-eigenfanction  i=l,2,...  of  A*A  such  that: 

A* A?  =  A J  n=l ,2, . . . 

n  n  n 


(see  [7]). 

A  sufficient  condition  on  D  for  ensuring  convergence  of  iteration  (27)  is 
the  following  (see  [21]): 


a . 

®^n  Pn^n ’ 

b. 

Pn  satisfies  0  <  pnkn  <  2 

for  all  n. 

c . 

Pn,  n=l,2,...  is  a  bounded 

sequence . 

It  is  interesting  to  remark  that  the  operator  AfcA  is  given  by  the  integral 
kernel:  sincp  and  therefore  PR,  n=l,2,...  are  the  prolate  spheroidal  wave 
functions  ( [9] ) . 


Many  operators  can  be  chosen  to  satisfy  conditions  a,  b  and  c.  In  [21], 
it  was  shown  that  it  is  sufficient  to  pick  D  =  G(AtA)  where  G(A)  is  a  polyno¬ 
mial  or  rational  function  such  that  0  <  AG(X)  <  2  for  0  <  1  (  1.  If  D  is  to  be 

2 

so  chosen,  (27)  will  converge  in  the  L  (2)  norm  to  the  solution  of  the  problem 
(Af)(x)  =  g(x),  x  6  F,  where  g:  Rn  ->  C  is  assumed  to  be  a  C-band-limited 
function.  If  we  now  apply  inverse  Fourier  transform  to  both  sides  of  (27)  we 
will  get  the  following  recursion 


gQ  =  0 

gn(x)  =  gn_1(x)  +  /q  e'",’WX  D{/ p  e  ^7Tiz’(g(z)  -  gn_^  (  z )  )  dz  )  ( w)  dw  (IS) 
Observe  that  when  D  =  I,  (28)  becomes  the  Landweber  Papou I i s - Gerchbe r g  algo¬ 


rithm.  Equation  (28)  shows  a  quite  general  version  of  this  classical 
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situation. 

In  the  remainder  of  this  section,  ve  will  present  some  numerical  simula¬ 
tion  results  comparing  the  generalization  (27)  with  the  classical  iteration 
(10').  To  this  end,  let  us  define 

7 

_  .  „  .  .  .sin  n/2  x. 

g:  R  ->  R:  g(x)  =  ( — — *  cos  nx 

The  Fourier  transform  of  this  signal  is  plotted  in  Fig.  1.  If  we  take  the 
interval  F  =  (-1,1)  as  the  known  part  of  g,  a  fairly  reasonable  reconstruction 
of  the  Fourier  transform  can  be  obtained  by  means  of  Discrete  Fourier 
Transform  (DFT)  of  129  samples.  This  result  is  plotted  in  Fig.  2.  It  is  clear 
that  the  two  peaks  are  easily  distinguished.  Or  <.„e  other  hand,  if  F  =  (-1/2, 
1/2)  the  situation  will  be  completely  different.  Figure  3  plots  the  result 
obtained  for  DFT  of  129  samples  in  (-1/2,  1/2).  This  means  that  restricting 
the  known  part  to  (-1/2,  1/2)  represents  an  irretrievable  loss  for  the  appli¬ 
cation  of  the  naive  inversion  technique.  In  other  words,  by  means  of  DFT  of 
samples  of  g  on  (-1/2,  1/2)  the  outstanding  features  of  the  spectrum  of  g  are 
lost.  Therefore,  we  think  that  g:  [-1/2,  1 / 2 J  ->  R  is  a  reasonable  test  exam¬ 
ple  for  our  numerical  simulations. 

We  first  apply  the  Landweber-Papoul i s-Gerchberg  iteration.  Figure  4a 
shows  the  very  poor  result  obtained  after  20  iterations.  Figure  4b  plots  the 
reconstructed  Fourier  transform  after  500  iterations.  In  this  case  the 

We  now  apply  the  more  general  procedure  given  by  (27)  for  three  different 

D’s. 

1.  D  =  (A* A  +  y I ) -1 
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For  this  operator  (27)  is  closely  related  to  the  Twomey-Tikhonov  regulariza¬ 
tion  method.  In  this  case,  iteration  (27)  becomes: 


fo  =  0 

(A*A  +  yl)fn  =  +  A*g,  nil  (29) 


It  is  worth  pointing  out  that  y  should  be 
yl  is  always  invertible. 

Figure  5a  shows  the  result  obtained 
and  Figure  5b  plots  the  reconstructed 
with  the  same  parameter  y.  In  both  cases 
ity. 


chosen  positive.  In  that  case,  A  A  + 

after  10  iterations  when  y  =  0.00005, 
Fourier  transform  after  20  iterations 
the  reconstructions  are  of  good  qual- 


Fizing  a  value  for  y  and  determining  the  number  of  iterat ; ons  are  by  no 
means  trivial  matters.  By  comparing  Figures  5a  and  5b  it  is  seen  that  the 
reconstruction  is  quite  sensitive  to  the  number  of  iterations  d  es.  We  think 
that  the  sensitivity  depends  also  on  the  parameter  y. 

Figure  6a  shows  the  result  after  10  iterations  obtained  by  applying  (29) 
when  y  =  0.005.  Figure  6b  plots  the  corresponding  result  for  y  =  0.005  and  50 
iterations.  By  comparing  Figures  5a  and  6a  it  is  seen  that  the  reconstruction 
is  very  sensitive  to  the  parameter  y  when  the  number  of  iterations  is  fixed. 

Figure  7  shows  the  reconstruction  obtained  for  y  =  0.000005  after  10 
iterations.  It  is  clear  that  for  a  fixed  number  of  iterations  the  smaller  the 
parameter  y  is,  the  more  distorted  (due  to  the  propagation  of  round-off 
errors)  the  reconstruction  will  be. 


In  spite  of  some  unanswered  questions,  the  main  conclusion  that  can  be 
drawn  trotn  these  examples  is  that  the  resolution  obtained  ir  Figures  5  i  .  1 
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6a  and  7  is  nncli  batter  than  that  of  Figares  4a  and  4b. 

2.  D  -  F(A*A)  and  FU)  -  804.375k6  -  3003X5  +  4504.5k4  -  3465k3  +  1443.75k2 
-  3151  +  31.5 

Some  reasons  for  choosing  such  a  D  are  veil  documented  in  [21].  We  think 
that  this  example  is  also  useful  to  realize  that  the  reconstruction  is  very 
sensitive  to  the  choice  of  D.  For  this  D,  Figure  8a  shows  the  result  after  10 
iterations.  It  is  seen  that  the  resolution  is  poor.  However,  Figure  8b  plots 
the  reconstructed  Fourier  transform  for  200  iterations  which  is  a  good  result. 
This  means  that  the  procedure  is  slower  compared  to  those  given  where  D  =  (A*A 

+  r I)"1. 

It  is  also  remarkable  that  by  using  a  fewer  number  of  iterations  than 
those  necessary  for  the  classical  Landweber-Papoul i s-Gerchberg  algorithm 
recursion  (27)  provides  a  better  reconstruction  (compare  Figures  4b  and  8b). 

3.  D  »  F(AtA)  and  FU)  -  l16. 

This  case  is  intended  to  be  an  example  where  the  speed  of  the  reconstruc¬ 
tion  seems  to  be  similar  to  that  of  the  classical  approach  (9)-(10'). 

Figure  9a  shows  the  reconstruction  obtained  after  500  iterations.  By  com¬ 
paring  9a  <  and  4b  is  is  seen  that  the  results  look  much  the  same.  Figure  9b 
plots  the  result  obtained  after  1,000  iterations.  By  comparing  9b  with  9a  it 
is  noticeable  that  the  reconstruction  of  the  Fourier  transform  was  improved  at 
the  cost  of  double  computational  effort. 

It  was  assumed,  so  far,  that  the  given  signal  is  not  contaminated  with 
any  noise.  Since  the  techniques  presented  in  1  and  2  above  represent  a  sub¬ 
stantial  improvement  of  the  classical  iteration  procedure  (8)-(10')  it  is 
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expected  that  the  noise  will  also  propagate  much  faster  in  the  reconstruction. 
Therefore,  a  stopping  rule  is  of  great  ioportance  for  practical  applications. 
It  is  also  important  to  analyze  what  is  the  performance  of  the  iteration  when 
the  known  range  of  g  is  smaller.  Some  related  examples  and  further  analyses 
are  given  in  [23]. 
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ABSTRACT 


This  paper  deals  with  the  numerical  continuation  problem  of  analytic  functions 

g(z)  given  by  g(z)  ■  /  K(z,t)  f(t)  dt  where  K  defines  a  bounded  operator  on 

2  ^ 

L  (Cl)  .  It  is  assumed  that  g  is  known  over  a  finite  segment  A  of  the  real  line 
where  g  is  to  be  sampled.  Our  continuation  techniques  emerged  partly  from 
a  generalization  of  the  Landweber  iteration.  We  show  that  a  certain  discrete 
approximation  of  the  proposed  iterative  technique  yields  two-step  (non¬ 
iterative)  algorithms  for  solving  the  continuation  problem.  We  also  prove 
convergence  of  these  approximations  to  the  sought  function  g.  Special 
emphasis  on  the  continuation  problem  for  the  case  K(z,t)  =  e  ~^2t  is  given 
and  some  related  numerical  examples  are  presented.  The  continuation  problem 
when  the  known  part  of  g(z)  is  contaminated  with  some  noise  is  addressed  and 
some  techniques  for  solving  this  problem  are  also  provided. 


I.  INTRODUCTION 


Lee  us  suppose  that  we  are  given  a  piece  of  a  n-dimensional  signal 
g(t),  t  ft  A  £Rn.  In  addition,  we  assume  that  g  is  obtained  from  some 
other  signal  f (x) ,  x  &Q  ,  through  a  linear  space  variant  system, 

,x)  f(x)  dx,  tSR",  (1) 

where  K  is  known.  The  goal  is  to  recover  the  ''real”  object  f(x),  x€:Q, 
from  a  finite  set  of  samples  of  the  "observed"  signal  g(t),  when  t  €:  A. 
This  problem  is  very  well  known  in  the  engineering  literature  ([1])  and  it 
has  been  extensively  studied  in  the  mathematical  literature  ([2]). 

A  very  important  case  is  obtained  from  (1)  when  K(t,x)  »  k(t-x), 
t,x  £]Rn.  In  that  case,  g  and  f  will  satisfy  the  following  relationship: 

g(u>)  -  k(u»)  •  Tf(oJ)  ,  Rn 


where  A  denotes  the  Fourier  transform  and  TF(x)  «  f(x)  if  x&Q,  and  0 
elsewhere.  It  Is  clear  that 


1  A  a 

1 -  “  TF (ui)  ,  for  all  u>  :  k(uj)  ^  0 


A 

Let  us  assume  that  k(uJ)  f  0  for  ui  6  N,  where  N  contains  a  non-void  open  set 

tl  a  A 

of  R  .  Since  TF  has  compact  support  (if  *1  is  compact)  then  TF  is  analytic; 

A  A 

so  the  knowledge  of  TF(uj),  when  uJ  £  N,  will  be  enough  to  determine  TF(uJ) 

for  any  other  uj£]Rn.  in  many  applications,  TF(**0,  tu£!Rn  will  describe 

by  itself  all  the  information  that  we  need  from  TF.  If  this  is  not  the 

A  n 

case,  then  we  should  proceed  to  compute  Tf  from  Tf(<o),  to  £  R  .  However, 


III.  2 


v  have  assumed,  so  far,  Chae  g(tu)  is  known  exactly  for  all  uJ6nn.  This 
will  not  be  the  case  if  g(t)  is  observed  on  the  set  A  1*11°,  or  on  a  finite 
subset  of  A  only.  This  shows  that  if  we  can  improve  our  knowledge  of  g 
(i.e.,  to  know  g(t),  when  t  A)  we  will  obtain  a  better  knowledge  of  g 
(i.e.,  to  compute  g(u>)  more  accurately). 

A 

In  many  cases,  k(uu)  ■  0  if  "  $ N,  where  N  is  assumed  to  be  compact. 
Therefore  g(«J)  -  0,  to  £  N,  which  assures  that  g  will  also  be  an  analytic 
function.  This  means  that  the  set  of  values  g(t),  t  €  A  will  determine 
g(t),  t  £  A.  This  shows  that  the  solution  f  to  the  equation 


g(z) 


/‘ 


k(z-x)  f(x)  dx 


can  be  approached  by  solving  a  continuation  problem  for  two  analytic  func- 

A 

tions:  g(z),  given  z  €  A,  and  TF(uj),  given  m6N.  It  is  important  to  notice 

A 

that  the  continuation  of  TF  can  be  stated  in  the  sense  of  equation  (1) 
since 


A 

TF  (o>)  ” 

and  e  plays  the  role  of  K(u/,x).  This  latter  continuation  problem 

shows  an  example  of  the  importance  of  including  space-variant  Kernels  in  our 
discussion. 

One  motivation  for  the  continuation  problem  we  give  above  is  the 
restoration  of  f.  However,  there  is  another  motivation:  In  many  cases,  we 
are  Interested  in  obtaining  knowledge  of  g(x) ,  when  x  £  A,  and  x  is  "close" 
to  A.  Some  examples  of  this  situation  are  known,  in  multidimensional  signal 


/• 


-2TTixu> 

e  f  (x) 


dx 


processing.  Let  us  suppose  we  apply  a  filter  to  a  given  image.  The  filter 
ideally  performs  over  an  infinite  extent  image .  However,  in  practice  we  are 
given  only  a  piece  of  the  image.  When  the  filter  is  applied  to  points  close 
to  or  on  the  border  of  the  real  image,  Inaccuracies  will  result,  if  we  assume 
some  arbitrary  numbers  for  the  unknown  values  of  the  image  (e.g.,  the  image 
is  assumed  to  be  periodic;  or  a  constant  number  is  assumed  for  the  unknown 
values).  The  filter  would  improve  its  performance  if  we  could  fill  out  the 
unknown  values  of  the  image  with  some  interpolated  information.  Thus, 
small  amounts  of  extrapolation  (i.e.,  to  extrapolate  a  small  region  beyond 
the  boundary  of  A)  can  be  of  great  help. 

In  the  following  sections  we  concern  ourselves  with  the  continuation 
of  the  function  g,  when  g  is  given  by  the  equation 

g(z)  -  f  K(z,t)  f(t)  dt,  z  €  (1) 

T; 

and  g(z)  is  known  for  z  €:  A  ClRn. 

One  possible  approach  for  getting  a  continuation  of  g  is  to  solve 
equation  (1)  in  terms  of  f  and  to  use  the  same  formula  (1)  to  obtain  the 
continuation  of  g.  We  will  use  the  Landweber-Strand  ([3], [4])  iterative 
procedure  to  get  the  solution  of  equation  (1) : 


fQ  *  initial  approximation. 


fn-.£n-l  +  DK  (S-KW  ' 


(’> 


where  K  denotes  the  adjoint  of  K  and  D  is  a  certain  operator  ([4]).  I: 
section  II,  we  will  extend  (2)  to  cases  where  K  is  not  compact  and  D  is 


suitable  positive  operator,  including  Strand's  iteration  ([4])  as  a  special 
case.  If  we  apply  K  to  both  sides  of  equation  (2)  we  will  get 

*n  "  *n-l  +  K  D  K*(s  -  Vl3  (3> 

gQ  “  initial  guess  , 

where  g  ■  Kf  is  a  function  defined  on  the  whole  . 
n  n 

Under  rather  general  conditions  for  K  and  f,  recursion  (2)  is  known  to 

be  convergent.  On  the  other  hand,  we  will  also  prove  that  the  sequence  gR 

approaches  g  uniformly  over  compact  sets  in  Cn.  Practical  computation  of  gn 

requires  equation  (3)  to  be  sampled;  this  means  that  discretization  of  g^f  g 
* 

and  KDK  are  unavoidable.  In  this  paper,  we  will  show  that  certain  natural 
discretization  provides  an  iterative  procedure  which  is  also  convergent  to  a 
sequence.  This  sequence  can  also  be  obtained  by  solving  a  system  of  equations 
where  D  plays  an  essential  role.  This  system  of  equations  provides  a  natural 
interpolation  for  the  sequence  by  means  of  an  analytic  function.  \'e  will 
also  prove  that  this  analytic  function  approaches  g  uniformly  over  compact 
sets  when  the  distance  between  samples  of  g  on  A  tends  to  zero. 

Thus,  we  will  obtain  a  reliable  technique  for  continuating  g(z),  because 
the  closer  the  samples  of  g  are  taken  on  A,  the  better  the  approximation 
to  g(z),  when  z  £A.  It  turns  out  that  our  continuation  technique  does 
not  require  the  computation  of  the  derivatives  of  g.  Section  IV  which  con¬ 
tains  the  main  results  is  derived  independently  of  Landweber-Strand ’ s  iteration 
However,  Landweber-Strand' s  iteration  is  the  origin  of  the  main  ideas  developed 
in  this  paper  and  therefore,  section  II  is  fully  devoted  to  this  iterative 
procedure  and  our  generalization.  Section  III  shows  the  relationships  between 


Landweber-Strand's  iteration  and  our  continuation  techniques,  which  are 

presented  in  Section  IV.  Section  V  includes  some  numerical  examples  for  the 
-2Wixv 

case  K(x,y)  -  e  .In  section  VI  we  present  the  continuation  problem 

for  the  case  in  which  the  observations  are  given  by  g(x)  -  g(x)  +  ^(x),  x  6  A 
Here,  denotes  a  continuous  perturbation  such  that  J^Cx)^  fe  for  some  £>0 
All  our  results  are  presented  for  the  one-dimensional  case.  However,  exten¬ 
sions  to  the  multidimensional  case  readily  follow. 

II.  LANDWEBER-STRAND’S  ITERATION 

In  this  section  we  present  Strand’s  iteration  in  a  more  general  setting 

2  2 

than  that  of  [4].  Let  K  :  L  (Cl)  •— L  (A)  be  a  bounded  linear  operator.  Let 

2  i 

g  €s  L  (A)  be  a  function  such  that  g  =*  +  g2»  where  g^  6  R(K)  and  g2  &  R(K) 

*  * 

(R(K)  denotes  the  range  of  K) .  In  addition,  let  D  :  R(K  )  -»R(K  )  be  a 

bounded  linear  symmetric  operator  which  is  assumed  to  be  one-to-one  (Vx  :  Dx 

* 

then  x  *»  0)  ,  and  such  that  DK  K  is  non-negative. 

Theorem  1 

Under  the  conditions  stated  above,  let  f  be  defined  as  follows 

n 

fa  ’  fn-!  +*DK‘<*  - 

2 

f  *  any  initial  guess  €  L  G  -)  (4) 

2  + 
Let  us  suppose  that  04  x  4  _____  .  Then,  f  converges  to  (I  -  P)f  +  f  , 

«  UK,  K.  ii  n  o 

where  f+  is  the  minimum  norm  least  squares  solution  of  Kf  =  g,  P  is  the 
orthogonal  projection  on  N(K)  =  j^v  £  L  (Q)  :  Kv  =  0  ]  and  I  denotes  the 
identity. 


This  theorem  is  a  simple  consequence  of  a  general  result  (Bialy  [5]) 
for  bounded  operators  on  Hilbert  spaces.  Let  us  recall  that  our  assumption 
g  £  R(K)  +  R(K)1"  ensures  that  the  problem  Kf  »  g  has  a  least  squares  solu¬ 
tion  ([6]).  It  is  also  known  that  the  set  of  least  squares  solutions  is 
(,x  :  K*Kx  *  K*g}.  This  means  that  we  need  to  pick  up  the  minimum  norm  solu¬ 
tion  for  the'  operator  equation 

K*Kx  »  K*g  (5) 


Since  D  is  supposed  to  be  one-to-one,  equation  (5)  is  equivalent  to 


DK*Kx  -  DK*g 


(6) 


Since  DK  K  *  H  is  a  symmetric  bounded  linear  and  non-negative  operator, 
Bialy' s  conditions  [5]  are  satisfied  and  therefore  the  iteration 


fn  *  Vl  +<(8'‘  “»-l> 

f  »  any  initial  approximation 


(7) 


2 

converges  in  the  L  (£1)  norm  to  (I  -  P)f 

< 

solution  for  the  consistent  equation  Hx 
It  is  also  clear  that  (7)  becomes 


+  f+,  where  f+  is  the  minim"m  norm 
1  g,  provided  that  0  4  *  4  j^^([5]) . 


f„  -  £„-l+  DK‘(s  - 


We  can  include  the  relaxation  parameter  into  D.  In  that  case,  convergence 
of  (4)  is  ensured  if  IlDK  fell  4  2.  Strand  ([4])  considers  equation  (4)  for 


some  special  D  and  assuming  K  to  be  a  completely  continuous  operator.  If 


the  latter  condition  is  satisfied  then  the  singular  value  theory  for  compact 
operators  will  apply  to  our  problem.  In  that  case,  if  <f>  denotes  the  family 
of  eigenfunctions  of  K  K,  for  n  :  1,  2,  and  ^  its  corresponding 

eigenvalues : 

K*KJ^  ■  ^n^n’  n  *  1*  2,  ... 

Strand's  conditions  on  D  become 

D<6  -  p  <6  ,  n  »  1,  2,  ...  (8a) 

n  n  n 

where  p  satisfies 
rn 

0  4  p  4  2  for  all  n,  (8b) 

*n  n 

p^,  n  =  1,  2,  ...  is  bounded  (8c) 

* 

It  turns  out  that  condition  (8a)  ensures  that  DK  K  is  non-negative.  It  is 

also  clear  that  (8b)  is  equivalent  to  ll DK  K  ft  4  2,  and  condition  (8c)  ensures 

that  D  is  a  bounded  operator  on  R(K  K) .  These  conditions  allow  for  a  partial 

study  relating  D  to  the  behavior  of  the  approximation  (4)  after  n  -iterations. 

Under  the  more  general  situation  stated  in  theorem  1,  some  rational 
* 

expressions  in  K  K  may  play  the  role  of  D.  Similar  analyses  to  those  of 
[4]  are  expected  for  this  situation.  However,  this  goes  beyond  the  scope 
of  this  paper.  In  what  follows,  we  will  restrict  our  attention  to  cases  in 
which  K  defines  a  bounded  integral  operator. 
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III.  APPLICATIONS  TO  THE  CONTINUATION  PROBLEM 


We  remarked  in  the  introduction  that  the  iterative  procedure  (2)  can 
be  used  for  getting  the  continuation  of  certain  analytic  functions  given  by 


the  formula 


g(z)  -  /  K(z,t)  f(t)  dt  ,  z6l 

Jn 

2  2 

where  K  defines  a  bounded  operator:  L  (ft)  —  L  (A)  and  g  is  observed  on  the 

set  A.  It  is  clear  that  if  we  define  g  »  Kf  ,  where  f  is  given  by  (2) , 

n  n  n 

we  will  obtain  the  recursion 


*n  '  *n-l  +  roK  Cg  '  «n-l> 

S0  ’  0  "  «„*  fo  -  0  « 

Since  f^  converges  to  f  where  f  is  the  minimum  norm  solution  of  the  equation 
g(z)  ■  /  K(z,t)  f(t)  ,  z  6  A  , 

JQ 

2  2 

and  since  the  convergence  is  in  L  (Cl) ,  then  g  converges  to  g  in  L  (A)  .  We 

n 

now  state  the  following  result: 


Corollars 


Let  us  suppose  that  K  satisfies  these  additional  properties: 

(1)  /  K(z,t)  h(t)  dt  defines  an  analytic  function  for  z  6  C  and  any 

h  6  L2(.U. 

f  2 

(2)  sup  /  |  K(z,  t)l  dt  4  C  <•  os  for  all  compact  sets  I  C  C. 


Then,  the  sequence  g  given  by  (9)  converges  uniformly  over  compact  sets  in 

n 

C  to  the  continuation  of  g(x),  x  6  A  for  complex  arguments  z  £  C. 

Proof 

Since  =  Kf^  is  also  an  analytic  function  over  the  complex  plane, 

2 

and  f^-*  f  L  (Q)  we  have  by  means  of  the  Schwartz  inequality: 

|(g  -  gn)(z)  U  ft  fn  -  fH2  *  [  f  l  K(z,  t)l  2  dt]^ 

for  €:  arbitrarily  small,  n  ^  n  (~)  and  P  a  compact  set  in  C. 

o 

Some  particular  operators  allow  a  bound  (as  the  type  of  that  given  by 
(2)  in  the  Corollary)  which  depends  only  on 

sup  ( Imz 1  . 
zGT 

In  that  case,  convergence  is  assured  to  be  uniform  over  the  real  line.  This 

is  the  situation  for  K(z,t)  =  e  Some  consequences  of  this  Corollary 

in  signal  processing  are  given  in  (7]. 

Realistic  applications  of  this  iterative  continuation  formula  involve 

★ 

discretization  of  g,  g^  and  the  operator  KDK  .  Applying  naive  quadratures 

it 

formulas  for  the  integral  K  (g  -  g^  we  will  obtain  the  following  recursion 
for  any  k  £  Z: 

N  ,  _  , 

S  (kA)  =  S  .  (kA)  I  KDK(iA,  • )  (kA)-(g(iA)  -  :  ,(iA))  (10) 

n  n-1  .  J  n-I 

i=-N 

K(:c,y)  is  the  complex  conjugate  of  K(x,y)  and  the  dot  in  K(iA,‘)  indicates 


IIlJLO 


the  variable  with  respect  to  which  D  is  being  applied  on  K;  &  is  a  positive 
number  used  to  sample  2N  +  1  times  the  function  g  on  (-a, a)  :  L  ■  and  e< 

is  a  relaxation  factor  independent  of  n.  It  would  be  desirable  that 
Sn(kfiO  *  gQ(ki) ;  however  it  is  obvious  that  the  sequence  SQ(klL),  keZ.  is 
not  obtained  by  sampling  gR.  Two  problems  arise  from  the  discrete  recursion 
(10).  The  first  is  related  to  the  convergence  of  the  sequence  when  n-»«  . 

The  second  problem  is  about  the  relationship  between  this  discrete  technique 
and  the  original  continuation  problem.  More  specifically,  the  quewt ion  is  whether 

lim  S  (k&) 
n 

n-»co 

will  approach  g  when  A—0.  Both  problems  are  addressed  in  the  next  section. 


IV.  ALGORITHMS  FOR  THE  CONTINUATION  PROBLEM 


In  what  follows,  we  will  assume  that  K(x,*)  6  L  (Cl)  for  all  x  £  A, 

*  2 
D  =*  D  is  positive  definite,  and  defined  over  the  whole  L  (A).  Under 

these  conditions  the  discrete  recursion  (10)  will  be  derived  by  means  of  a 

rather  different  approach.  Let  us  consider  the  following  matrix  L 


L  :  Jlij ,  -N  4  i,  j  L  N. 

l±j  -  &(KDK(iA,-)  )(jA)  (1 

It  is  simple  to  prove  that  L  is  a  hermitian  matrix,  i.e.  1.,  ■  i. .  for  all 

a-J 

i,j.  If  we  assume  that  K  has  the  following  additional  property: 


Y  c  K(ii,x)  =  0,  for  all  x  fe  Cl  -  c .  —  0  ,  i  =  -N,  . .  .  ,  NT 
i— N  1 


€ 
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Chen  L  will  be  positive  definite  (without  condition  (12)  L  is  non-negative) . 
This  shows  that  the  system  of  equations 

Lc  -  d  , 

will  always  have  a  unique  solution  for  any  d  €  C^+^.  In  particular,  we  can 
find  i  ■  -N,  N  that  satisfy 

N  _ 

A  V  (KD  K(jA,*)  )  (iA)  Y,  -  (Kf)(i A),  -N  4  i  4  N  (13) 

j»-N  2 

It  is  also  worth  pointing  out  that  after  Y.  is  found  in  (13)  we  will  have  a 

3 

natural  interpolation  formula  for  g(iA),  by  means  of  an  analytic  function, 
i.e. , 

N  __ 

g  (z)  »  A  Y  (KD  K(jA,*)  )(z)Y  (14) 

A  jhN  J 

Let  us  show  how  this  approach  relates  to  recursion  (10) .  Since  L  is  a  posi¬ 
tive  definite  matrix,  then  ^ ,  -N  4  j  4  N  can  be  obtained  by  means  of  the 
following  well-known  iterative  procedure: 

Y^n)  =  Y^n-1)  +  o((g(iA)  -  L/n_1))  (15) 

o(  a  certain  positive  relaxation  parameter  and 

Y^  =  lim  Y^n^,  for  all  i  £  [-N,N]. 
n— oo 

We  now  apply  A (KD  K(iA,*)  )  (kA)  ,  k €  Z,  ,  to  both  sides  of  formula  (15)  to  get 


N 


iA  y  (KD  K(iA,  •)  )(kA).Y,(n)  =  A  T  (KD  K(iA,  * )  )  (kA)  . 


(n-1) 


i*-M 


i=-N 


+  °c  A  2  (KD  K(iA,*)  )  (kA)  .  (g(iA)  -  LY(n_1)) 


(16) 


III. 12 


We  have  shown,  so  far,  that  sampling  iteration  (g)  produces  a  point- 
wise  convergent  sequence  which  can  be  computed  by  an  exact  technique.  This 
procedure  consists  of  two  steps.  We  first  solve  a  system  of  linear  equations 

(13),  obtaining  i  -  -N . N;  then,  we  use  to  get  an  analytic 

function  g^(z) ,  z  6  C  which  interpolates  the  set  of  samples  g(iA),  i  - 
where  iA  6 A  for  all  i  6  [-N,N].  This  function  is  given  by  equation  (14): 

N 

g.(z)  -  A  V-  (KD  K(j&,*)  )(z).Y 

A  j-N  j 

The  remainder  of  this  section  is  devoted  to  showing  the  relationship  between 
g^(z)  and  g(z).  We  will  need  some  lemmas. 

Lemma  1 

Let  (hm)m  be  a  family  of  analytic  functions  over  Cn  and  let  us  assume 

that  (h  )  is  uniformlv  bounded  over  compact  subsets  of  Cn  (i.e., 
m  m 


for  all 
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a 

i 


■ 


compact  sacs  P £  Cn  chare  exists  a  constant  c  :  sup  Ih  (z) 1  4  c  for  all  m) . 

IE1  '  r 


Under  these  conditions  there  exists  a  sequence  h  ,  k  — oo  and  certain  analytic 

“k 

function  h  such  that  h  converges  to  h  uniformly  over  compact  subsets  of  Cn. 

®k 


For  a  proof  of  this  well-known  result  see,  for  instance,  [8].  In  what  fol- 

<  2a 

lows,  we  will  show  that  g  ,  -  ■  2^+1  satisfy  the  hypothesis  of  lemma  1, 
when  D  and  f  satisfy  some  additional  conditions.  We  will  also  need  another 
previous  result  from  the  optimization  theory  in  Hilbert  spaces. 

Lemma  2 

Let  (H,  <d,>  )  be  a  Hilbert  space.  Let  x^,  -N  4  i  4  N  be  2N+1  linear 
independent  elements  of  H  and  c^,  -U  4  i  4  N  be  2N+1  complex  numbers.  Then, 
the  following  optimization  problem: 


minimize  4x,x> 

subject  to  x  G  H 

4x,xi>  =  c^,  -N  i  4,  N 


is  uniquely  solved  by 


o 


x 


N 


x. 

1 


(17a) 


(17b) 


(13) 


where  ft  are  determined  by  solving  the  system  of  equations: 


L  -  cj  •  -sil‘s 


(19) 


For  a  proof  of  this  lemma  see,  for  instance,  [9].  We  are  new  able  tc  prove 


Ill  .14 


Theorem  2 

Let  us  suppose  that  D  satisfies:  there  exists  m  >  0  :  <£x,Dx>  >,  m^t,x> 

2  2  a 

for  all  x  £  L  ,  and  assume  that  there  exists  a  function  t  €.  L  such  that  DC  ■  f. 

2a 

In  that  case,  the  family  of  analytic  functions  g^,  given  by  formulas 

(13)  and  (14)  is  uniformly  bounded  over  compact  subsets  in  C. 


Proof 


will  use  lemma  2  where  x^  -  A**  K(iA, •)  and  c^  -  g(iA),  -N  4  14  N. 


The  underlying  Hilbert  space  is  L  (0)  equipped  with 


<x,y>  -  /  x(uj)  Dy(u>)  d<u  . 


Under  these  conditions,  the  optimization  problem  (17a)  becomes 


minimize  .  /  h(ui)  (Dh)  (u))  dui 
2 

subject  to  h  £  L  (Q) 


A*5  •  J h(u>)D^K(ii,  *)j  («»)  du)  »  g(ik) 


for  -N  4  1  4  N 


(20a) 


(20b) 


It  turns  out  that  the  solution  to  that  problem  is 


*°(w)  -  K*  t  A  K(iA,uJ)  , 

i— N  1 


where  ^  satisfies 


AT  MKOA.cu),  K(iA,«i)> 

*  zr*_  \r  J  s 


g(iA)  ,  -N  4  i  4  N 
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By  dafinition  of  <  ,  £,  condition  (22)  becomes 


N 

A,  T  P.  (KD  K(jA,‘)  )(iA)  -  g(iA)  ,  -N  4  1  u 

1— N  J 


This  last  system  of  equations  coincides  with  (13)  and  therefore,  ft  »  ^ , 

-N  ^  j  4  N.  It  is  also  important  to  notice  that  g^can  be  written  as  follows 


g(z)  -  aV*0(uJ)  D  K(z ,  • )  (uO  &u> 

L  JC1 


(23) 


where  is  the  optimum  given  by  equation  (21).  Since  d£  *  f,  for  certain 
tfe  L^,  the  function  A*1.  €  is  a  feasible  point  for  the  optimization  prob- 
lem,  that  is  to  say,  A  .Z  satisfies  condition  (20b).  Hence, 


*°(o>)  Dl° (ul)  duj  4  A'1 1  €  dcu 


(24) 


'Cl 


because  K  is  the  "minimun-norm"  function  among  those  which  satisfy  (20b) . 
Now,  the  sought  property  of  the  family  g^  readily  follows  from  (23)  and  (24) 
From  (23)  we  obtain 


UA(z)U  J*.11*0II2.Hd  k(z ,  • ) 


(25) 


Using  equation  (24) ,  and  the  assumptions  made  on  D  and  K,  (25)  becomes 

|g^(z)l  4  (fZ(v)  dZ(uj)  d*)*  HDll2.(/  1  K(z ,  t)  |  2  dt)'5  .  (25') 

If  we  now  assume  that  z£P  a  compact  set  in  C,  equation  (25')  becomes 


1  g^(z)l  4  c-p 

where  c^  depends  only  on  T. 


for  all  z  €.  T1  , 
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We  ere  now  eble  to  epply  Lemma  1  to  obtain 


Theorem  3 


g.(z) 

A 


-  &  f  *,<KD  KCJA.O  >(*>  *  *  6  c  .  A-  -2^. 

j— N  J 


be  the  analytic  function  obtained  with  the  condition 


-N  4  1  4  N  :  A  £  *  (KD  K(jA,‘)  )  (iA)  -  g(iA)  -  (Kf )  (iA)  (13) 


j— N 


Let  us  assume  that  D  and  K  satisfy  all  the  conditions  stated  In  Theorem  2. 

In  that  case,  there  exists  a  subsequence  A  —  0  where  m  —  <x>  such  that  g. 

m  n 

m 

approaches  g  uniformly  over  compact  sets  inC. 


Proof 


By  lemma  1,  there  exists  a » subsequence  A  -»  0  and  an  analytic  function 


h(z),  such  that  g  approaches  h  uniformly  over  compact  sets  inC.  In 

ul 

addition,  h  *  g  because  g  is  a  family  of  equicontinuous  functions  (a  simple 

consequence  of  the  uniform  boundness)  and  therefore,  since  g„  (iA )  •  g(iA  ) 

a  m  m 

when  -[2a/A  ]  A.  i  4  [2a/A  ]  and  the  set  of  numbers  iA  :  0  4  i4<a, 
mm  m 

-[2a/Am]  414  [2a/Am]  is  dense  in  (-a, a)  we  conclude  h(z)  *  g(z)  for 
z  (-a, a).  Since  both  functions  are  analytic,  the  same  identity  holds  for 


We  would  like  to  remark  on  some  important  points  and  to  give  some 


examples  before  ending  this  section. 
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II 

« 


The  property  that  DK  K  is  positive  was  not  used  in  the  proof  of  Theorem  3. 
On  the  contrary,  this  property  is  required  for  the  recursive  computation  of  f 
by  means  of  formula  (2).  It  is  worth  pointing  out  that  we  have  used,  so  far, 
equally  spaced  samples  of  g(iA)  where  A  -  2^+1  *  However,  similar  results 
to  those  of  theorem  2  and  theorem  3,  and  different  formulas  for  interpolating 
samples  (as  (14))  can  be  obtained  for  non-equally  spaced  samples  and  for 
other  regular  distributions.  It  is  also  easy  to  see  that  we  could  have  used 
the  following  sampling  sets: 

1  iV  >  0ik4“!  -\4 14  V  wtt  } 

k  k  k 


ft 


where  is  a  non-decreasing  infinite  sequence  of  natural  numbers.  In  that 

case,  by  means  of  theorem  3  we  would  obtain  a  subsequence  of  N^,  k  ■  1,  2,  ...» 

say  ©  ,  n  =*  1,  2,  ...,  such  that  g^  approaches  g  uniformly  over  compact  sets 
n  % 

in  the  complex  plane.  If  we  repeat  this  procedure  for  every  non-decreasing 

infinite  sequence  of  natural  numbers  we  will  conclude  the  convergence  of 

2a 

g.  to  g,  where  A„  =  -  ,  and  N  is  any  natural  number.  This  observation 

0^  N  2N+1  - — 

is  important  because  it  means  that  our  approximation  to  g  over  a  compact  set, 

by  means  of  gA  will  be  good  for  all  N  > N  if  it  is  good  for  N  .  This  obser- 
N  00 

vation  also  applies  to  more  general  sequences  of  sampling  sets. 

We  now  give  some  specific  examples  of  operators  which  may  play  the  role 
of  D  in  the  continuation  technique.  Probably  the  simplest  case  is  D  =  I,  the 
identity  operator.  In  that  case,  the  interpolation  formula  (14)  becomes 


g  ,(2) 
A 


•  f—  \ 

-  y  K^K(j A,-);  (z)  Y. 

j=-N  J 


» 

* 


1 


(26) 
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where  are  chosen  to  satisfy  g^(iA)  ■  g(iA),  1  «  -N,  . ..,  N.  Another 
*  ,  -1 

example  is  D  ■  (K  R  +  <I)  ,  where  o<  is  a  real  positive  number.  More  gen¬ 

erally,  if  P,^  denotes  the  .polynomial 

k  . 

P.  (X)  “  T.  a<  *  »  where  >  0  ,  a.  >  0 

*  i-0  1  1  1 

,  * 

then  D  -  Pfe(K  K)  will  satisfy  all  the  conditions  required  by  theorem  3  and 

-1  2 
so  does  D  .  The  examples  given  above  satisfy  R(D)  -  L  (ft) .  This  property 

2 

ensures  that  f  £  R(D)  for  any  f  6  L  (ft)  whatsoever.  In  particular  if  f  is 
such  that  Kf  ■  g,  we  will  have  f  &  R(D).  Hence,  no  relationships  between  f 
and  D  are  needed  to  apply  theorem  3.  Unfortunately,  this  is  not  always  the 
case.  The  other  property  which  D  is  required  to  satisfy  is 


3  m  >  0  ;  4Dx,x  >  >  m<x,x>  for  all  x  6  L  (ft) 


(27) 


This  assumption  is  by  no  means  mild.  Simple  positive  operators  may  not  have 
this  property.  For  example,  let  us  define 


D  :  L2(Q)  -  L2(n)  :  (D  ;.)  (w)  -  h(<u).£(co) 


(28) 


where  h  is  a  bounded,  real  and  positive  function.  Even  if  f  6  R(D) ,  that  is 
to  say 

]£! 


a 


(28') 


D  may  not  satisfy  (27)  and  therefore,  the  hypotheses  of  theorem  3  may  not  be 
satisfied.  However,  the  thesis  of  that  theorem  still  holds.  To  show  that, 
let  us  assume 

*  l2 


'Q 


TT 

h 


(29) 
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which  is,  in  fact,  a  weaker  condition  than  (28).  Then,  the  analytic  function 
g  can  be  written  as 


(“ft  "  8  .  ^l2(0) 

h  h 

where  D*4  :  L2(fl)  —  L2(0)  :  (T)\)  (ou)  *  tfyaJ)  .£(«*>) .  We  now  call 
and  we  apply  theorem  3  when  K  is  replaced  by  H  (notice  that  H  is 
operator)  and  when  no  D  is  present  (i.e.,  D  *  I).  We  obtain  the 


(26): 


N  _ 

g .  (z)  -  A  V  H  H(jA,0  (z)  * 

*  J— N 


H  -  KD  “ 

an  integral 

interpolation 


or  equivalently 


g  (z) 


l 


N  _ 

A  V  KD  K(jA,  •)  (z) 
j=-N  J 


(30) 


which  Is  the  formula  we  should  obtain  in  theorem  3  applied  to  this  situation. 

The  great  flexibility  in  choosing  D  shown  in  theorem  3  may  be  useful 
in  two  different  ways: 

1.  Some  more  a  priori  information  about  g  (or  f)  may  suggest  picking 
up  some  D's  to  get  closer  interpolation  formulas.  The  numerical 
evidence  of  the  next  section  gives  an  example  about  this  situation. 


2.  The  noise  problem.  When  g  is  contaminated  with  some  noise,  the 
effect  of  D  may  be  of  great  importance,  especially  when  the  system 
of  equations  (13)  is  solved  iteratively.  This  has  been  shown  to  be 
the  case  when  we  deal  with  the  iterative  procedure  (2)  for  computing 
f  (see  [I]).  Some  numerical  properties  of  related  techniques  for 
solving  the  problem  in  the  presence  of  noise  were  recently  reported 
in  [13]  when  K(x,t)  =  e 


To  finish  this  section  we  would 


like  to  ooint  out  that  another 
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interesting  question  remains  unanswered:  What  are  the  additional  conditions 
that  £  and  D  should  satisfy  to  get  an  approximation  result  such  as  that  of 
theorem  3  when  D  is  any  positive-definite  linear  operator? 


IV.  NUMERICAL  EXAMPLES  AND  DISCUSSION 

We  will  consider  in  this  section  some  numerical  examples  for  the  case 
K(z,t)  «  e“^Tt*2t.  The  special  interest  in  this  particular  operator  is  well 
known  and  we  tried  to  explain  some  of  the  reasons  in  Section  I. 

In  this  situation,  interpolation  formula  (14)  becomes 


(31) 


The  case  D  *  I  was  addressed  in  [12].  However,  in  ref.  [12]  it  was  also 

* 

assumed  that  g  €:  R(K  K)  to  derive  the  corresponding  interpolation  formula 
and  no  relationships  between  g^and  g  were  shown.  On  the  contrary,  we  do  not 
need  any  further  restriction  on  g  than  g  €  R(K) ,  that  is  to  say. 


g(z) 


/„■ 


f (t)  e"2Uizt  dt  ,  z  6  A 


and  since  the  kernel  satisfies  all  the  properties  required  by  theorem  3  we 

also  know  that  g  g  uniformly  over  compact  sets.  In  ref.  [10]  it  was 

& 

shown  that  for  this  particular  situation  the  convergence  of  g ,  to  g  is  uniform 

A 

over  the  whole  real  line. 

r  if  i2 

If  we  assume  that  D  is  given  by  (28)  and  if  /  da)  <  co  where 

A  a  _  XI  h 

h  :  1R  -»]R  is  bounded  and  h(oi)  «  0  if  oj  £  then  we  can  apply  formula  (30) 


to  obtain: 
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N 

8  (z)  -AT  *  h(z-jA)  (32) 

A  j—  N  J 

In  this  situation,  the  system  of  equations  we  need  to  solve  to  compute  <S\  , 
j  =  -N,  N  has  the  following  form 

N 

g(iA)  =  A  £  tf  h((i-j)A)  (33) 

jj— ~N  J 

It  turns  out  that  the  matrix  involved  in  (33)  is  Toeplitz  and  therefore, 
2 

0((2N  +  1)  )  algorithms  can  be  used  for  solving  the  system.  We  have  developed 
an  iterative  approach  for  solving  (33)  whose  relaxation  parameter  does  not 
depend  on  the  number  of  points  inside  A  (for  A.  fixed),  when  A  is  varying 
([11]). 


Our  numerical  experience  will  be  based  on  the  following  functions: 


g1(z) 


(— cos  1  z 


g2(z) 


sin  2TT(z-^) 
TT(z-Js)  ' 


sin  2TT(z+4) 


z  £  C 


g3(z) 


/  sin  TT  z  \ 

V  tr2  /  * 


z  6  C 


It  is  clear  that  these  functions  are  analytic  and  they  can  be  written  as 


with  f  6  L"  (Cl) 


f(x) 


-2tTizx  . 
e  dx 


JQ 

and  Q=  (-1,  1). 


The  interval  (-0.5,  0.5)  is  assured  as 


:ne 


.1122 


set  of  known  values.  We  will  take  33  equally  spaced  samples  In  (-0.5,  0.5) 
with  A  ”  ~  ,  for  each  case.  The  continuation  will  be  given  by  means  of 
65  equally  spaced  samples  of  g^(z)  taken  In  (-1,  1),  corresponding  to  the 
points  z  ■  iA,  1  £^[-32,  32]. 

In  the  first  case,  g^(z)  Is  given  by1  the  integral 

f  0,5  -  ,  .  -2lTio)z  . 

I  f.  (to)  e  dto 

2-0.5  1 


where  f  is  drawn  in  figure  l.a.  The  set  of  33  samples  of  g^^  is  given  in  table 
l.a.  Figure  l.b  shows  the  function  gx  in  (-1,1).  These  samples  provide  a  very 
poor  description  of  the  spectrum  of  g(x),  x6E,  because  their  discrete  Fourier 
transform  (DFT)  does  not  distinguish  the  corresponding  pair  of  peaks.  The  DFT 
is  shown  in  table  l.b  and  in  figure  2.2.  Compare  this  with  the  DFT  obtained 
when  the  33  samples  are  taken  in  (-1,1)  (fig.  2.b).  In  this  case  we  have 
chosen  D  =>  I  and  therefore,  g  is  given  by  (32)  with 


h(x) 


sin  21fx 
TTx 


Table  l.c  shows  the  values  of  the  continuation  g^"^  (iA)  ,  32  > li|  ^  17, 

A  ■  ■—  .  Table  l.d  shows  the  real  values  g^(iA),  32  >  [i  \  ^  17.  It  is 

seen  that  the  continuation  is  of  good  quality.  We  now  compute  the  DFT  of 

the  sequence  defined  by  the  given  samples  g^CiA),  0  4  I  i  l  ^  16  plus  the 

continuation  samples  g^(iA),  32  >  [il  17.  The  result  is  given  in  table 

A 

1.3.  Notice  that  the  two  peaks  are  now  distinguished  very  clearly.  This 
example  shows  the  importance  and  effectiveness  of  the  continuation  procedure. 

The  second  case  is  intended  to  be  an  example  of  the  effectiveness  of 
the  algorithm  when  the  main  peaks  of  the  function  to  be  continuated  are  not 
completely  included  into  the  known  range  (fig.  3).  In  this  example,  we  use 
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0 

a 


* 

D  ■  K  K  +  -'I,  where  »  0.5.  The  function  obtained  by  means  of  the  formula 
KD(K(x, •)) (z)  can  be  written  in  this  case  as 


where 


sinc1>0(z-y) 


sinc^g(y-x)  dy  +  0.5  sin^  g(z-x) 


sinc1.0(uj) 


sin  2TTuj 
TTuJ 


(34) 


D 


t 


The  function  obtained  by  the  continuation  procedure  is  sampled,  giving 
(2) 

g^'(i&),  32  >lii  >  17  .(table  2.b).  The  values  of  the  known  samples  of 
g^  are  listed  in  table  2. a,  and  the  real  values  ,  32  ^  [  i  1  ^  17  are 

given  in  table  2.c.  The  conclusion  that  can  be  drawn  by  comparing  tables  2.b 
and  2.c  is  that  the  continuation  procedure  performs  a  fairly  reasonable 
extrapolation.  It  is  worth  noticing  the  accuracy  obtained  at  the  main  peaks. 
However,  these  peaks  are  still  close  to  the  known  range  of  the  given  function. 
Further  numerical  evidence  is  necessary  when  the  location  of  the  peaks  is 
more  critical. 

The  third  case  is  to  provide  a  rough  comparison  of  the  performance  of 
two  different  D's.  The  function  (fig.  4)  is  to  be  continuated  from  its  samples 

g^CiiO,  0  4  l  i  l  4  16  by  means  of  *  sinc^  g(x-y)  and  D,  given  by  (34). 

Table  3. a  shows  the  known  samples.  Tables  3.b  and  3.c  list  the  values  of 
both  continuations.  Table  3.d  contains  the  real  values  of  g^  taken  on  the 
same  points  at  which  the  continuations  are  sampled.  It  is  seen  that  the 
continuation  provided  by  is  closer  to  g^  than  that  of  on  most  of  the 
sample  points.  This  phenomenon  is  not  yet  fully  explained. 

Probably  the  main  conclusion  that  may  be  drawn  from  these  examples  is 


> 


that  33  samples  of  the  given  function  was  always  enough  to  get  a  quite  reason' 
able  continuation.  It  goes  without  saying  that  this  result  Is  due  to  some 
peculiar  properties  of  the  given  functions  g^»  gr}  and  g^»  in  terms  of  the 
chosen  D’s.  A  deeper  study  of  these  properties  might  shed  more  light 
on  the  role  played  by  D  in  the  continuation  technique. 

IV.  CONTINUATION  AND  NOISE 

We  have  assumed,  so  far,  that  our  samples  of  g,  taken  on  A,  are  not 
contaminated  with  any  noise.  However,  a  more  realistic  model  should  consider 
that  the  given  observations  are  g(x),  x  6  A,  where  g  “  g  +  *|.  It  is  clear 
that  g  and  r]  cannot  be  separated.  In  general,  is  not  a  smooth  noise  and 
therefore,  we  need  to  change  the  formulation  of  the  continuation  problem 
because  g  is  not  a  piece  of  some  analytic  function.  We  will  assume  that 
rj  :  A  — C  is  a  continuous  function  which  satisfies  [rj(x)l  4  €  ,  for  all  x  €  A, 
where  €  is  some  known  constant.  Under  these  assumptions,  we  could  seek  for 
a  function  h(z)  which  is  given  by  f K(z,t)  s(t)  dt,  for  all  z  6  C,  s 
and  such  that  ] g(x)  -  h(x)(  4  £ ,  for  all  x  €.  A.  It  is  obvious  that  this 
problem  has  a  solution.  However  it  may  not  be  unique  and  therefore,  h  will 
be  different  from  g.  This  means  that  In  the  absence  of  any  other  a  priori 
information  about  g,  the  continuation  of  g  from  g  becomes  a  de^cate  matter. 

Among  the  infinitely  many  functions  which  may  play  the  role  of  h 
(i.e.,  which  satisfy  the  two  conditions  stated  above),  we  may  put  some  addi¬ 
tional  constraint  to  the  problem  to  guarantee  a  unique  solution:  for  all 

€  >  0.  The  only  restriction  we  should  impose  to  the  additional  constraint 
(in  the  absence  of  some  more  _a  priori  information  about  g)  is  that  the 


solution  approaches  (in  some  convenient  way)  the  true  function  g  when 
tends  to  zero.  This  convergence  must  involve  the  values  of  g  oeyond  A. 

In  this  section,  we  develop  a  technique  which  will  use  the  following  three 
constraints: 


Vz) 


IS{ 


s^-(t)  K(t,z)  dt  ,  2  €:  C 


|  h  (x)  -  g(x)  |  A  6  ,  for  all  x  6  A 


s_  -z  L  minimizes  <Cs,s>_ 


fS 

Jn 


s (t) (Ds) (t)  dt 


(35a) 

(35b) 

(35c) 


among  those  which  satisfy  (35a)  and  (35b) 


In  formula  (35c),  D  denotes  a  linear,  bounded  operator  which  satisfies 
<Dx,x^  m<x,x>  for  all  x,  with  m  >  0.  We  will  give  a  procedure  to  compute 
h^  and  we  will  show  that  h^  approaches  g  uniformly  over  compact  sets  in  the 
complex  plane  when  g(z)  =  /k(z,  t)  f(t)  dt,  z  €  C  and  d£  *  f  for  some  t  €  L2(  }) 


Q 


Theorem  4 


Let  K  be  an  integral  operator  satisfying  the  properties  stated  in 
theorems  2  and  3.  Let  D,  g  and  g  be  as  defined  above,  and  let  us  assume  that 
£  =  D?  for  a  certain  2.  €  L  (D)  .  If  we  take  A^  =  ,  for  all  N  ^  1,  and  if 

we  define 


N 


h  ,  (z)  -  An  £  Y  .  KD  K(j^,-)  (z)  ,  z  6 


(36) 


j=-N 


where  Y^ *  j  =  -S , 


. ,  S  is  the  optimum  of  the  following  problem: 


minimize 


Z  ■  K 


j— N  J 


(37a) 


subject  to: 
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Hi  KD  KU^.O  ay  -  .  -8^  i<  H 

I  -  g(i^)  U  €  .  -N  6  1  4  N 


(37b) 


(37c) 


Then,  there  exists  a  subsequence 


approaches  hg  uniformly  over 


compact  sets  in  t£ie  complex  plane,  where  hg  is  defined  by  conditions  (35a), 
(35b)  and  (35c). 


Proof 


We  will  follow  the  same  ideas  as  those  of  theorem  3.  It  is  clear  that 

h  is  well-defined  because  condition  (35b)  defines  a  closed  convex  subset 
2 

of  L  (CD  •  It  is  also  clear  that  h.  is  well-defined  and  that  condition  (37c) 

% 

can  be  rewritten  as 


-NiiiN  :  1  h 


A  •  8(1V  I  4  e 


It  is  also  possible  to  rewrite 


Z  Vj 

j— N  J  J 


as  follows: 


JLvj  ■  3t  <■« 

■  /  J/j  K0V>  •  k?(kV->) (t)  dt 

T> 

N  N 

•<  z  ^K(jA,o  ,  £  r  K(jA .')> 

j—  N  J  N  j—  N  2  N  D 
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and  therefore,  the  optimization  problem  can  be  restated  by  using  lemma  2 
as  follows 


minimize  <s,s>_ 


(38a) 


subject  to 


s  ~L  : 


-N  4  i  4  N  ,  l<s,  L*  K (i\.*)>D  -  g(iAN)  |  4  £  (38b) 


It  is  also  important  to  notice  that  h,  can  be  wirtten  in  the  following  way: 


hA  (z)  -  A^/sj(t)  •  d[k(z,-)]  (t)  dt 

N  Jr,  ‘ 


(39) 


where  is  the  optimum  solution  to  the  problem  (38a) ,  (38b)  which  is  given 
by  the  formula 


SN(t)  "  AN  I  VC)  ’  C  6  Q 

j«-N  J 


By  means  of  (39) ,  we  obtain 


lh^(z>[  4  a*h  •  Hsjll2  .  [Id  k(z,-)  || 


(40) 


n  n  2  -h 

Since  f  ■  Dt,  c  €:  L  (Q)  ,  then  A?2  also  satisfies  condition  (38b). 
(40)  and  the  other  properties  assumed  for  D  and  K  we  conclude 

:z)  U  [  IU t)  d£(o  dzj 

"N 

where 


:rom 


|  hL  (z)  |  4  (  i( t)  De(t)  dt)  .1  Dll2  .  -t  .  C  , 
u 


C_  =  sup 
*  z£ 


P  (  f[K(z.OS2  dt)'2  , 

JQ 


a  compact  set  in  C. 


We  now  apply  lemma  1  to  the  sequence  h,  ,  >i  0  to  get  an  analytic  functic 


> 


g~  which  is  a  uniform  limit  on  compact  sets  of  a  certain  subsequence  hA  > 

6  ^  ** 
k  >  0.  It  is  clear  that  gfe  satisfies  j  g^(x)  -  g(x)(  4  €  for  all  x  €  A 

2a  - 

because  h^  is  an  equic on tinuous  family,  [  iA^,  -N  4  i  4  N»  is  a 

dense  subset  of  A  and  is  continuous.  Therefore  g^  satisfies  (35b). 

Since 

(z)  -  A^yS^(t).D^K(z,*>]  (t)  dt. 


and 


a 


ilA^  Sj#2  4  c,  N>1 


then 


gc(z) 


■/« 

•h 


(DS?)(t)  K(z.t)  dt  , 


2  W 

where  S°  is  the  weak  limit  in  L  (Si)  of  a  certain  subsequence  of  AlT  .S“  . 

\  k 

This  means  that  gg  satisfies  equation  (35a) . 

h  m 

Since  is  a  weak  limit  of6^mS*m  ,  k  ^  0,  a  subsequence  of  Nfc, 


then 


I 


S“(t)  DS“(t)  dt  4  sup 
€  €  k 

n  a 


,/4»k 


(t)  DSjn(t)  dt 


(41) 


In  addition,  by  definition  ofS®, 

N 


kN  fSU 


(t)  DS“(t)  dt  4 


6  fsj 

Jn~ 


(t)  DS^(t)  dt  ,  for  all  N  >  0. 


and  by  (41)  we  obtain 


I3 


'Is 


S*(t)  DS*(t)  dt  4  j  S£(t)  DSfe(t)  dt 


(42) 


By  (42)  and  (35c)  we  conclude 


Since  satisfies  condition  (35b)  then  we  have 

Si  -  S_  a.e.  in  Q.  • 

Hence,  h£(z)  *  gfe(z)  for  all  z  and  therefore  h^  is  given  as  the  limit  of 

h^  ,  uniformly  over  compact  sets  in  C. 

Nk 

It  can  be  proved  that  there  exists  a  sequence  of  positive  numbers 

£  0  such  that  h£  —  g  uniformly  over  compact  sets  in  the  complex  plane. 

n 

The  proof  can  be  done  by  means  of  similar  arguments  to  those  of  theorem  3  and 
theorem  4  and  it  will  not  be  included  here. 

To  finish  this  paper,  we  would  like  to  point  out  that  efficient  algo¬ 
rithms  for  solving  the  optimization  problem  (37a)  -  (37b)  -  (37c)  are  being 
studied.  We  hope  to  present  these  numerical  results  in  a  future  paper. 
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